17 最終課題：KADAI1・KADAI2

7. コピペ用Pythonスクリプト

以下のコードは、配布する KADAI01_chl_kmeans_fill.py と KADAI02_soi_nino_persistence_fill.py の中身をそのまま貼り付けたものである。学生は必要な方をコピーして、Jupyter Lab の新しいPythonファイルまたはNotebookセルに貼り付けて実行する。

7.1 KADAI01：Chl季節変化パターン分類

# =========================================================
# KADAI01: 10年分Chlデータから季節変化パターンをK-meansで分類する
# 穴埋め版スクリプト
# =========================================================

from pathlib import Path
import numpy as np
import pandas as pd
import xarray as xr
import matplotlib.pyplot as plt

from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler

# ---------------------------------------------------------
# 0. 設定
# ---------------------------------------------------------
DATA_FILE = "GlobColour_CHL1_monthly_100km_AV_201601_202512.nc"
VAR_NAME = "CHL1"
K = 5

OUT_DIR = Path("output_KADAI01")
FIG_DIR = Path("figures_KADAI01")
OUT_DIR.mkdir(exist_ok=True)
FIG_DIR.mkdir(exist_ok=True)

# 太平洋域の切り出し範囲
LON_MIN, LON_MAX = 120, 290
LAT_MIN, LAT_MAX = -60, 60

# ---------------------------------------------------------
# 1. NetCDFを読み込む
# ---------------------------------------------------------
ds = xr.open_dataset(____)   # ヒント: DATA_FILE
print(ds)

# Chl変数を取り出す
chl = ds[____]               # ヒント: VAR_NAME

# 0以下の値はlog10が取れないので欠損値にする
chl = chl.where(chl ____ 0)   # ヒント: >

# log10(CHL)に変換する
logchl = np.____(chl)         # ヒント: log10

# ---------------------------------------------------------
# 2. 10年分から12か月気候値を作る
# ---------------------------------------------------------
# time.monthごとに平均する
logchl_clim = logchl.groupby(____).mean(____, skipna=True)
# ヒント1: "time.month"
# ヒント2: "time"

print(logchl_clim)

# ---------------------------------------------------------
# 3. 太平洋域を切り出す
# ---------------------------------------------------------
# lonが0〜360の場合を想定
sub = logchl_clim.sel(
    lon=slice(____, ____),    # ヒント: LON_MIN, LON_MAX
    lat=slice(____, ____)     # ヒント: LAT_MIN, LAT_MAX
)

# 緯度が降順の場合に備える
if sub.sizes.get("lat", 0) == 0:
    sub = logchl_clim.sel(
        lon=slice(LON_MIN, LON_MAX),
        lat=slice(LAT_MAX, LAT_MIN)
    )

print(sub)

# ---------------------------------------------------------
# 4. K-means用の特徴量行列を作る
# ---------------------------------------------------------
# shape: month, lat, lon -> lat, lon, month
arr = sub.transpose("lat", "lon", "month").values
nlat, nlon, nmon = arr.shape

# 各格子点を1行、12か月を12列にする
X = arr.reshape(____, ____)   # ヒント: nlat*nlon, nmon

# 欠損がある格子点を除く
valid = np.all(np.isfinite(X), axis=____)  # ヒント: 1
X_valid = X[valid, :]

print("Number of valid grid cells:", X_valid.shape[0])

# ---------------------------------------------------------
# 5. 標準化してK-means
# ---------------------------------------------------------
scaler = StandardScaler()
X_scaled = scaler.fit_transform(____)     # ヒント: X_valid

kmeans = KMeans(n_clusters=____, random_state=0, n_init=20)  # ヒント: K
labels_valid = kmeans.fit_predict(____)   # ヒント: X_scaled

# 元の格子に戻す
labels_flat = np.full(nlat * nlon, np.nan)
labels_flat[valid] = labels_valid + 1  # クラスタ番号を1始まりにする
labels_map = labels_flat.reshape(nlat, nlon)

# ---------------------------------------------------------
# 6. クラスタマップを描く
# ---------------------------------------------------------
plt.figure(figsize=(11, 5))
pc = plt.pcolormesh(sub["lon"], sub["lat"], labels_map, shading="auto", cmap="tab10")
plt.colorbar(pc, label="Cluster number")
plt.xlabel("Longitude")
plt.ylabel("Latitude")
plt.title("K-means clusters of seasonal log10(CHL) patterns")
plt.grid(True, alpha=0.3)
plt.savefig(FIG_DIR / "KADAI01_cluster_map.png", dpi=300, bbox_inches="tight")
plt.show()

# ---------------------------------------------------------
# 7. 各クラスタの平均季節変化を描く
# ---------------------------------------------------------
months = np.arange(1, 13)
summary_rows = []

plt.figure(figsize=(9, 5))
for c in range(1, K + 1):
    idx = labels_valid == c
    mean_cycle = np.nanmean(X_valid[idx, :], axis=0)
    plt.plot(months, mean_cycle, marker="o", label=f"Cluster {c}")

    max_month = int(months[np.nanargmax(mean_cycle)])
    min_month = int(months[np.nanargmin(mean_cycle)])
    amplitude = float(np.nanmax(mean_cycle) - np.nanmin(mean_cycle))

    summary_rows.append({
        "cluster": c,
        "n_grid_cells": int(idx.sum()),
        "max_month": max_month,
        "min_month": min_month,
        "seasonal_amplitude_log10chl": amplitude,
        "mean_log10chl": float(np.nanmean(mean_cycle))
    })

plt.xlabel("Month")
plt.ylabel("Mean log10(CHL)")
plt.title("Mean seasonal cycle for each cluster")
plt.xticks(months)
plt.grid(True, alpha=0.3)
plt.legend()
plt.savefig(FIG_DIR / "KADAI01_cluster_mean_seasonal_cycle.png", dpi=300, bbox_inches="tight")
plt.show()

# ---------------------------------------------------------
# 8. 結果をCSVとNetCDFで保存する
# ---------------------------------------------------------
summary = pd.DataFrame(summary_rows)
summary.to_csv(OUT_DIR / "KADAI01_cluster_summary.csv", index=False)
print(summary)

cluster_da = xr.DataArray(
    labels_map,
    coords={"lat": sub["lat"], "lon": sub["lon"]},
    dims=("lat", "lon"),
    name="cluster"
)
cluster_da.to_netcdf(OUT_DIR / "KADAI01_cluster_map.nc")

print("Done: KADAI01")

7.2 KADAI02：SOIとNiño3.4の持続性比較

# =========================================================
# KADAI02: SOIとNiño3.4の自己相関から持続性を比較する
# 穴埋め版スクリプト
# =========================================================

from pathlib import Path
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

SOI_FILE = "SOI_timeseries_updated.txt"
NINO_FILE = "sstoi.indices"

OUT_DIR = Path("output_KADAI02")
FIG_DIR = Path("figures_KADAI02")
OUT_DIR.mkdir(exist_ok=True)
FIG_DIR.mkdir(exist_ok=True)

START_DATE = "1982-01-01"
END_DATE = "1997-12-01"
MAX_LAG = 36

# ---------------------------------------------------------
# 1. 関数
# ---------------------------------------------------------
def standardize(x):
    """平均0、標準偏差1に標準化する。"""
    x = np.asarray(x, dtype=float)
    return (x - np.nanmean(x)) / np.nanstd(x)


def corr_no_index_alignment(a, b):
    """indexの自動整列を避けて相関係数を計算する。"""
    a = np.asarray(a, dtype=float)
    b = np.asarray(b, dtype=float)
    mask = np.isfinite(a) & np.isfinite(b)
    if mask.sum() < 3:
        return np.nan
    return np.corrcoef(a[mask], b[mask])[0, 1]


def autocorrelation(x, max_lag=36):
    """lag 0〜max_lagまでの自己相関を返す。"""
    x = np.asarray(x, dtype=float)
    out = []
    for lag in range(max_lag + 1):
        if lag == 0:
            r = ____          # ヒント: lag 0の自己相関は1.0
        else:
            r = corr_no_index_alignment(x[:-lag], x[____:])  # ヒント: lag
        out.append(r)
    return np.array(out)


def first_lag_below(acf, threshold):
    """自己相関がthreshold以下になる最初のlagを返す。"""
    for lag in range(1, len(acf)):
        if acf[lag] <= ____:  # ヒント: threshold
            return lag
    return np.nan

# ---------------------------------------------------------
# 2. SOIを読む
# ---------------------------------------------------------
soi = pd.read_csv(SOI_FILE, sep=r"\s+")
print(soi.head())
print(soi.columns)

# 年月からdateを作る。列名はデータに合わせて確認すること。
soi["date"] = pd.to_datetime(dict(
    year=soi[____],     # ヒント: "YEAR"
    month=soi[____],    # ヒント: "MONTH" または "MO"
    day=1
))

soi_df = soi[["date", ____]].rename(columns={____: "SOI"})
# ヒント: SOI列名を指定する

# ---------------------------------------------------------
# 3. Niño3.4を読む
# ---------------------------------------------------------
# sstoi.indices は空白区切りを想定する。
# 列構造が違う場合は print(nino.head()) と print(nino.columns) で確認する。
nino = pd.read_csv(NINO_FILE, sep=r"\s+", comment="#")
print(nino.head())
print(nino.columns)

# ここでは YEAR, MONTH, NINO34 という列名に整える想定。
# 必要に応じて列名を修正する。
nino["date"] = pd.to_datetime(dict(
    year=nino[____],    # ヒント: 年の列
    month=nino[____],   # ヒント: 月の列
    day=1
))

nino_df = nino[["date", ____]].rename(columns={____: "NINO34"})

# ---------------------------------------------------------
# 4. 共通期間にそろえる
# ---------------------------------------------------------
merged = pd.merge(____, ____, on=____, how=____)
# ヒント: soi_df, nino_df, "date", "inner"

merged = merged[(merged["date"] >= START_DATE) & (merged["date"] <= END_DATE)]
merged = merged.dropna(subset=[____, ____])
# ヒント: "SOI", "NINO34"

print(merged.head())
print(merged.tail())
print("Number of months:", len(merged))

# ---------------------------------------------------------
# 5. 標準化
# ---------------------------------------------------------
merged["SOI_std"] = standardize(merged[____])       # ヒント: "SOI"
merged["NINO34_std"] = standardize(merged[____])    # ヒント: "NINO34"

# 共通時系列を保存
merged.to_csv(OUT_DIR / "KADAI02_common_standardized_timeseries.csv", index=False)

# ---------------------------------------------------------
# 6. 時系列図
# ---------------------------------------------------------
plt.figure(figsize=(12, 4))
plt.plot(merged["date"], merged["SOI_std"], label="SOI standardized")
plt.plot(merged["date"], merged["NINO34_std"], label="Niño3.4 standardized")
plt.axhline(0, linewidth=0.8)
plt.xlabel("Year")
plt.ylabel("Standardized anomaly")
plt.title("Standardized SOI and Niño3.4 time series")
plt.grid(True, alpha=0.3)
plt.legend()
plt.savefig(FIG_DIR / "KADAI02_standardized_timeseries.png", dpi=300, bbox_inches="tight")
plt.show()

# ---------------------------------------------------------
# 7. 自己相関
# ---------------------------------------------------------
lags = np.arange(0, MAX_LAG + 1)
acf_soi = autocorrelation(merged["SOI_std"].values, max_lag=MAX_LAG)
acf_nino = autocorrelation(merged["NINO34_std"].values, max_lag=MAX_LAG)

autocorr_df = pd.DataFrame({
    "lag_months": lags,
    "SOI_autocorrelation": acf_soi,
    "NINO34_autocorrelation": acf_nino
})
autocorr_df.to_csv(OUT_DIR / "KADAI02_autocorrelation.csv", index=False)

plt.figure(figsize=(10, 5))
plt.plot(lags, acf_soi, marker="o", label="SOI")
plt.plot(lags, acf_nino, marker="o", label="Niño3.4")
plt.axhline(0, linewidth=0.8)
plt.axhline(0.5, linestyle="--", linewidth=0.8, label="r = 0.5")
plt.xlabel("Lag (months)")
plt.ylabel("Autocorrelation")
plt.title("Autocorrelation of SOI and Niño3.4")
plt.grid(True, alpha=0.3)
plt.legend()
plt.savefig(FIG_DIR / "KADAI02_autocorrelation.png", dpi=300, bbox_inches="tight")
plt.show()

# ---------------------------------------------------------
# 8. 持続性指標
# ---------------------------------------------------------
summary = pd.DataFrame({
    "index": ["SOI", "NINO34"],
    "lag_r_below_0p5": [first_lag_below(acf_soi, 0.5), first_lag_below(acf_nino, 0.5)],
    "lag_r_below_1overE": [first_lag_below(acf_soi, 1/np.e), first_lag_below(acf_nino, 1/np.e)],
    "lag_r_below_0": [first_lag_below(acf_soi, 0.0), first_lag_below(acf_nino, 0.0)],
})

summary.to_csv(OUT_DIR / "KADAI02_persistence_summary.csv", index=False)
print(summary)

# ---------------------------------------------------------
# 9. 発展：相互相関（余裕がある人のみ）
# ---------------------------------------------------------
# 正のlag: SOIがNiño3.4に先行する、と解釈する。
# 必須ではないが、余裕があれば考察に使ってよい。

print("Done: KADAI02")