これまで学んだ操作を、実データを使って確認する。ファイルを読む、中身を見る、時刻を作る、共通期間をそろえる、関数を書く、図とCSVを保存する、という流れを最終課題につなげる。
この講義では、Pythonの文法そのものだけでなく、海洋・気候データを実際に扱うための一連の流れを学んできた。最終課題では、その流れを自分で組み立てる。
このページでは、説明を読むだけでなく、短いコードを実行しながら復習する。最終課題のコードは長く見えるが、実際にはここで確認する基本操作を組み合わせたものである。
データ解析では、いきなり計算に入らない。まずファイルを読み、データの形を確認し、解析できる形に整える。
まず、解析で使うライブラリを読み込む。図やCSVを保存するためのフォルダも作っておく。
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from pathlib import Path
import xarray as xr
out_dir = Path("review_outputs")
out_dir.mkdir(exist_ok=True)
print("準備完了")まずは表形式の時系列データを読む。ここでは SOI_timeseries_updated.txt を使う。
soi = pd.read_csv("SOI_timeseries_updated.txt", sep=r"\s+")
print(soi.head())
print(soi.columns)
print(soi.shape)列名を確認したら、SOIの時系列図を描く。列名が自分のファイルと違う場合は、下の YEAR, MONTH, SOI を実際の列名に合わせる。
soi["date"] = pd.to_datetime(dict(
year=soi["YEAR"],
month=soi["MONTH"],
day=1
))
soi = soi.dropna(subset=["SOI"])
plt.figure(figsize=(12, 4))
plt.plot(soi["date"], soi["SOI"])
plt.axhline(0, color="k", linewidth=0.8)
plt.xlabel("Year")
plt.ylabel("SOI")
plt.title("SOI time series")
plt.grid(True, alpha=0.3)
plt.savefig(out_dir / "review_soi_timeseries.png", dpi=300, bbox_inches="tight")
plt.show()2つの時系列を比較するときは、同じ年月どうしを比べる必要がある。ここではSOIとNiño3.4を日付で結合し、共通期間だけを残す。
Niño3.4のファイル sstoi.indices は、列名や区切りがデータによって違う可能性がある。まず中身を確認する。
nino = pd.read_csv("sstoi.indices", sep=r"\s+", comment="#")
print(nino.head())
print(nino.columns)
print(nino.shape)列名を確認した上で、日付列を作る。ここでは例として、年の列を YR、月の列を MON、Niño3.4の列を NINO34 としている。
# 必要に応じて、自分のファイルの列名に合わせて変更する
nino["date"] = pd.to_datetime(dict(
year=nino["YR"],
month=nino["MON"],
day=1
))
nino = nino[["date", "NINO34"]].dropna()日付をキーにして、SOIとNiño3.4を結合する。
merged = pd.merge(
soi[["date", "SOI"]],
nino[["date", "NINO34"]],
on="date",
how="inner"
)
merged = merged.dropna(subset=["SOI", "NINO34"])
print(merged.head())
print(merged.tail())
print(merged.shape)SOIとNiño3.4は単位も値の範囲も違う。そのまま重ねると、どちらの変動が大きいのか比較しにくい。そこで標準化を行う。
def standardize(x):
x = np.asarray(x, dtype=float)
return (x - np.nanmean(x)) / np.nanstd(x)
merged["SOI_std"] = standardize(merged["SOI"])
merged["NINO34_std"] = standardize(merged["NINO34"])
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, color="k", linewidth=0.8)
plt.xlabel("Year")
plt.ylabel("Standardized anomaly")
plt.title("Standardized SOI and Niño3.4")
plt.legend()
plt.grid(True, alpha=0.3)
plt.savefig(out_dir / "review_soi_nino34_standardized.png", dpi=300, bbox_inches="tight")
plt.show()同じ処理を何度も使う場合は、関数にしておくとよい。関数を作ると、コードが短くなり、間違いも見つけやすくなる。
自己相関では、同じ時系列を時間方向にずらして相関を計算する。この処理はSOIにもNiño3.4にも使うので、関数にする。
def corr_no_index_alignment(a, b):
"""pandasの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]自己相関は、同じ時系列を時間方向にずらして比較する方法である。ラグが大きくなっても自己相関が高い場合、その時系列は過去の状態を長く保っていると考えられる。
Pythonでは、配列の番号は1からではなく0から始まる。要素が n 個ある場合、最初の要素は x[0]、最後の要素は x[n-1] である。
x = [10, 20, 30, 40, 50]
x[0] = 10
x[1] = 20
x[2] = 30
x[3] = 40
x[4] = 50ここで x は、SOIやNiño3.4の値を時間順に並べた配列である。
x = [x[0], x[1], x[2], ..., x[n-1]]x[:-1] は最後の1個を除いた配列、x[1:] は最初の1個を除いた配列である。
x[:-1] = [x[0], x[1], x[2], ..., x[n-2]]
x[1:] = [x[1], x[2], x[3], ..., x[n-1]]したがって、lag = 1 では次のように「ある月の値」と「1か月後の値」を比べている。
x[0] と x[1]
x[1] と x[2]
x[2] と x[3]
...
x[n-2] と x[n-1]def autocorrelation(x, max_lag=36):
x = np.asarray(x, dtype=float)
r_list = []
for lag in range(max_lag + 1):
if lag == 0:
r = 1.0
else:
r = corr_no_index_alignment(x[:-lag], x[lag:])
r_list.append(r)
return np.array(r_list)max_lag = 36
lags = np.arange(max_lag + 1)
soi_ac = autocorrelation(merged["SOI_std"], max_lag=max_lag)
nino_ac = autocorrelation(merged["NINO34_std"], max_lag=max_lag)
acf_df = pd.DataFrame({
"lag_month": lags,
"SOI_autocorr": soi_ac,
"NINO34_autocorr": nino_ac
})
acf_df.to_csv(out_dir / "review_autocorrelation.csv", index=False)
print(acf_df.head())plt.figure(figsize=(10, 5))
plt.plot(lags, soi_ac, marker="o", label="SOI")
plt.plot(lags, nino_ac, marker="o", label="Niño3.4")
plt.axhline(0, color="k", linewidth=0.8)
plt.axhline(0.5, linestyle="--", label="r = 0.5")
plt.xlabel("Lag (months)")
plt.ylabel("Autocorrelation")
plt.title("Autocorrelation of SOI and Niño3.4")
plt.legend()
plt.grid(True, alpha=0.3)
plt.savefig(out_dir / "review_autocorrelation.png", dpi=300, bbox_inches="tight")
plt.show()KADAI1では、10年分の月別Chlデータから、12か月の気候値を自分で作る。ここではNetCDFを読み、変数と次元を確認する。
ds = xr.open_dataset("GlobColour_CHL1_monthly_100km_AV_201601_202512.nc")
print(ds)変数名を確認したら、Chlを取り出す。ここでは変数名を CHL1 とする。
chl = ds["CHL1"]
print(chl.sizes)
print(chl.time.values[0], chl.time.values[-1])Chlは0以下にならないため、0以下の値は欠損として扱う。その後、log10変換する。
chl = chl.where(chl > 0)
logchl = np.log10(chl)10年分の月別データから、12か月気候値を作る。
chl_clim = logchl.groupby("time.month").mean("time", skipna=True)
print(chl_clim)
print(chl_clim.sizes)jan = chl_clim.sel(month=1)
plt.figure(figsize=(10, 5))
plt.pcolormesh(jan["lon"], jan["lat"], jan, shading="auto")
plt.colorbar(label="log10(CHL)")
plt.xlabel("Longitude")
plt.ylabel("Latitude")
plt.title("January climatology of log10(CHL)")
plt.savefig(out_dir / "review_chl_january_climatology.png", dpi=300, bbox_inches="tight")
plt.show()K-meansでは、各データ点を特徴量ベクトルとして扱う。KADAI1では、各格子点の12か月気候値を特徴量にする。
つまり、K-meansは各格子点の緯度・経度そのものではなく、12か月のChl季節変化パターンが似ているかどうかで分類する。
X = np.column_stack([
chl_clim.sel(month=m).values.ravel()
for m in range(1, 13)
])
print(X.shape)解析結果は、Jupyter上で見るだけでなく、ファイルとして保存する。
acf_df.to_csv(out_dir / "autocorrelation.csv", index=False)plt.savefig(out_dir / "figure.png", dpi=300, bbox_inches="tight")cluster_da.to_netcdf(out_dir / "kmeans_cluster_map.nc")10年分の月別Chlデータを読み、12か月気候値を自分で作成し、K-meansで海域を分類する。
NetCDFxarray月別気候値log10変換K-meansPNG/CSV出力SOIとNiño3.4を共通期間にそろえ、標準化し、自己相関を計算する。どちらがより長い持続性をもつか、またその差は大きいと言えるかを考察する。
CSV/TXTpandas共通期間標準化自己相関PNG/CSV出力| 確認項目 | 内容 |
|---|---|
| 入力データ | 指定されたファイルを読み込んでいるか。 |
| 中身の確認 | head, columns, shape, print(ds) などで確認したか。 |
| 前処理 | 欠損値、log変換、標準化、共通期間抽出などを適切に行っているか。 |
| 解析 | K-meansまたは自己相関の処理が正しく実行されているか。 |
| 図 | タイトル、軸ラベル、凡例、カラーバーがあるか。 |
| 出力 | PNGとCSVが保存されているか。 |
| 考察 | 図や表から読み取れることを、自分の言葉で説明しているか。 |