Pythonデータ解析のおさらい

1. これまで何をやってきたか

この授業では、Pythonの文法そのものを暗記することよりも、実際の海洋・気候データを使って、データ解析の流れを身につけることを重視してきた。

最初は、タヒチとダーウィンの海面気圧から南方振動指数（SOI）を計算した。その後、SOIの移動平均、自己相関、フーリエ変換、スペクトル解析、SSTとの相関、フィルタ処理、ENSOコンポジット、さらに海面高度偏差（SLA）から地衡流や渦度を計算するところまで進んだ。

2. Python解析の共通パターン

どの回でも、ほぼ同じ順番で解析している。大事なのは、いきなり計算式を書くのではなく、まずデータを読み、形を確認し、欠損値を処理し、必要な部分を取り出してから解析することである。

典型的なコードの流れ

# 1. ライブラリを読み込む
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# 2. データを読む
df = pd.read_csv("data.csv")

# 3. 中身を確認する
print(df.head())
print(df.columns)
print(df.shape)

# 4. 欠損値を処理する
df = df.dropna()

# 5. 必要な列を取り出す
x = df["TIME"].to_numpy()
y = df["SOI"].to_numpy()

# 6. 解析する
y_mean = np.mean(y)
y_anom = y - y_mean

# 7. 図にする
plt.plot(x, y_anom)
plt.xlabel("Time")
plt.ylabel("Anomaly")
plt.show()

この流れは、CSVでもNetCDFでもほとんど同じである。ただし、CSVでは pandas、NetCDFでは xarray を使うことが多い。

3. CSVデータのおさらい

SOI、NINO index、時系列データの多くは、CSVやテキスト形式の表データとして扱った。このようなデータでは pandas を使う。

CSV・テキストデータの読み込み

import pandas as pd

df = pd.read_csv("SOI_timeseries_updated.txt", sep=r"\s+")

pd.read_csv() はCSVだけでなく、空白区切りやタブ区切りのテキストデータも読むことができる。空白区切りのときは sep=r"\s+" を使った。

DataFrameとは何か

pandasで読み込んだ表は DataFrame と呼ばれる。Excelの表に近いものだと思えばよい。行と列があり、列名を使ってデータを取り出せる。

print(df.head())      # 最初の5行を表示
print(df.columns)   # 列名を表示
print(df.shape)     # 行数・列数を表示

列を取り出す

time = df["TIME"]
soi  = df["SOI"]

df["SOI"] は、DataFrameの中からSOIという列だけを取り出す操作である。

欠損値を取り除く

df = df.dropna(subset=["SOI"]).reset_index(drop=True)

欠損値があると、平均、相関、回帰などの計算結果が nan になることがある。そのため、計算前に欠損値を取り除く。

期間を選ぶ

df2 = df[(df["YEAR"] >= 1951) & (df["YEAR"] <= 1997)]

角括弧 [] の中に条件を書くと、その条件を満たす行だけを取り出せる。この操作は、解析対象期間をそろえるときに非常に重要である。

pandasからNumPy配列へ変換する

x = df["SOI"].to_numpy()

pandasの列はそのままでも計算できるが、自己相関やFFTなどでは NumPy 配列に変換して扱うことが多い。

4. NetCDFデータのおさらい

SSTやSLAのような地図データは、CSVのような2次元の表ではなく、時間・緯度・経度を持つ多次元データである。このようなデータは、NetCDF形式で配布されることが多い。

NetCDFの読み込み

import xarray as xr

ds = xr.open_dataset("sla_data.nc")
print(ds)

xr.open_dataset() でNetCDFファイルを開くと、Dataset が得られる。Datasetには、緯度、経度、時間、SLA、ADT、流速など、複数の変数が入っている。

変数を取り出す

sla = ds["sla"]
lat = ds["latitude"]
lon = ds["longitude"]

CSVでは df["SOI"] のように列を取り出した。NetCDFでは ds["sla"] のように変数を取り出す。

範囲を切り出す

sub = ds.sel(
    longitude=slice(135, 150),
    latitude=slice(25, 40)
)

sel() は座標値を使ってデータを選ぶ操作である。たとえば、黒潮周辺だけを切り出したい場合に使う。

条件でマスクする

sla2 = sla.where(np.isfinite(sla))

where() は条件を満たす場所だけを残し、条件を満たさない場所をNaNにする。海だけ残す、欠損を除く、特定のしきい値以上だけ残す、という処理でよく使う。

CSVとNetCDFの違い

5. 時系列解析のおさらい

SOIやNINO3.4のように、時間とともに変化するデータを時系列という。時系列解析では、値そのものだけでなく、周期性、遅れ、平滑化、周波数成分を見る。

移動平均

移動平均は、前後の値を平均して、短い時間スケールのギザギザを弱める方法である。

df["SOI_13mo"] = df["SOI"].rolling(window=13, center=True).mean()

偏差と標準化

気候データでは、値そのものよりも「平年からどれだけずれているか」を見ることが多い。これを偏差、または anomaly という。

clim = df["SOI"].mean()
anom = df["SOI"] - clim

さらに標準偏差で割ると、単位の違うデータを比較しやすくなる。これを標準化という。

z = (df["SOI"] - df["SOI"].mean()) / df["SOI"].std()

自己相関

自己相関は、同じ時系列を少しずらして、自分自身とどれくらい似ているかを見る方法である。周期性がある時系列では、特定のラグで相関が高くなる。

x1 = soi[:-lag]
x2 = soi[lag:]

num = np.sum(x1 * x2)
den = np.sqrt(np.sum(x1**2) * np.sum(x2**2))
r = num / den

lag は時間のずれを表す。月データなら lag=12 は12か月のずれである。

フーリエ変換とスペクトル

時系列を時間方向に見るだけでは、何か月周期の変動が強いのかは分かりにくい。フーリエ変換を使うと、時系列を周期成分に分解できる。

x = soi - np.mean(soi)
X = np.fft.fft(x)
power = np.abs(X)**2

6. 関係を見る解析

SOIとNINO3.4、SOIとSSTのように、2つのデータの関係を調べるときには、相関係数、回帰直線、決定係数、ラグ相関などを使う。

相関係数

相関係数は、2つの変数がどれくらい一緒に増減するかを表す数値である。値は -1 から 1 の範囲をとる。

mask = np.isfinite(x) & np.isfinite(y)
r = np.corrcoef(x[mask], y[mask])[0, 1]

mask は、xとyの両方が欠損していない場所だけを選ぶための条件である。

回帰直線

回帰直線は、xが変化したときにyがどのように変化するかを、直線で近似する方法である。

a, b = np.polyfit(x[mask], y[mask], 1)
y_fit = a * x + b

決定係数 R²

決定係数は、回帰直線がyの変動をどれくらい説明できるかを表す。

y_obs = y[mask]
y_pred = a * x[mask] + b

ss_res = np.sum((y_obs - y_pred)**2)
ss_tot = np.sum((y_obs - np.mean(y_obs))**2)
r2 = 1 - ss_res / ss_tot

相関係数が高いほど、回帰直線で説明できる割合も大きくなる。ただし、相関が高いことは、必ずしも因果関係を意味しない。

ラグ相関

ラグ相関は、2つの時系列の片方を時間方向にずらしながら相関を計算する方法である。どちらの変化が先に起こるのかを考える手がかりになる。

lags = np.arange(-24, 25)
rs = []

for lag in lags:
    if lag < 0:
        x1 = x[:lag]
        y1 = y[-lag:]
    elif lag > 0:
        x1 = x[lag:]
        y1 = y[:-lag]
    else:
        x1 = x
        y1 = y

    mask = np.isfinite(x1) & np.isfinite(y1)
    rs.append(np.corrcoef(x1[mask], y1[mask])[0, 1])

7. 空間データ解析のおさらい

SSTやSLAのようなデータは、時間だけでなく、緯度・経度方向にも広がっている。このようなデータでは、地図として描くだけでなく、空間的な勾配や流速を計算することができる。

2次元マップを描く

plt.pcolormesh(lon, lat, sla, shading="auto")
plt.colorbar(label="SLA (m)")
plt.xlabel("Longitude")
plt.ylabel("Latitude")
plt.show()

pcolormesh は、緯度経度格子上の値を色で表示する関数である。SST、SLA、渦度、EKEなどの地図表示で使う。

緯度経度を距離に直す

緯度経度は角度であり、メートルではない。空間微分を計算するには、1度が何メートルに相当するかを考える必要がある。

Re = 6371000.0
meters_per_deg_lat = np.pi * Re / 180.0
meters_per_deg_lon = meters_per_deg_lat * np.cos(np.deg2rad(lat))

経度1度の距離は緯度によって変わる。赤道付近では長く、高緯度では短くなる。

SLAの勾配から地衡流を求める

海面高度に傾きがあると、地衡流が生じる。SLAを η とすると、地衡流は次のように表される。

相対渦度とEKE

流速が求まると、流れの回転の強さや、渦に伴う運動エネルギーも計算できる。

相対渦度は、流れが時計回り・反時計回りにどれくらい回転しているかを表す。EKEは、流速が強い場所で大きくなり、黒潮や渦の活動が強い場所を示す指標になる。

ベクトルを重ねる

plt.pcolormesh(lon, lat, sla, shading="auto")
plt.quiver(lon[::3], lat[::3], u[::3, ::3], v[::3, ::3])
plt.colorbar(label="SLA (m)")
plt.show()

quiver を使うと、地図上に流速ベクトルを重ねることができる。ただし、すべての格子点に矢印を描くと見づらくなるので、[::3] のように間引くことが多い。

8. 最後に：どの解析でも同じ考え方

ここまで扱ったテーマは、SOI、ENSO、SST、SLA、地衡流などさまざまである。しかし、Pythonでデータ解析を行うときの基本的な考え方は共通している。

最低限覚えておきたいこと

• CSVは pandas、NetCDFは xarray で読む。
• 解析前には、必ずデータの中身、列名、次元、欠損値を確認する。
• 平均を引くと偏差、標準偏差で割ると標準化になる。
• 相関は「似ているか」、回帰は「どれくらい説明できるか」を見る。
• 時系列では、移動平均、自己相関、ラグ、スペクトルが重要である。
• 空間データでは、緯度経度、距離、勾配、地図表示を意識する。
• 最後は必ず、図が何を意味しているのかを自分の言葉で説明する。