09 海面水温の長期トレンド

1. 背景：全球海面水温は一様に変化するのか

全球平均の海面水温は長期的には上昇傾向を示すが、海のどこでも同じ速さで上がっているわけではない。ある海域では強く上昇し、別の海域では変化が小さかったり、期間によっては低下して見えることもある。

そこで今回は、全球海面水温データを格子点ごとに時間方向へ線形回帰し、海面水温のトレンドを地図として可視化する。また、求めた傾きが統計的に 0 と区別できるかどうかを p 値で判定する。

全球海面水温トレンド（1982–2024）の例。赤は昇温、青は冷却を示す。
出典： Copernicus Marine Service, GLOBAL_OMI_TEMPSAL_sst_trend

Copernicus 図の Classification

2. 使用データ

本演習では NOAA Extended Reconstructed SST version 5（ERSSTv5）の月平均海面水温データを使う。ファイルはこのページと同じディレクトリに置いてあるものを使う。

ERSSTv5 の仕様

3. 計算の流れ

1. ERSSTv5 の NetCDF ファイルを読み込む
2. 解析したい期間を切り出す
3. 時間を「年の小数」に変換する
4. 各グリッドで y = ax + b の線形回帰を行う
5. 傾き a をトレンドとして保存する
6. 傾きの p 値を求め、有意な格子点を判定する
7. トレンドマップと有意性マップを描く

4. Pythonでよく出る機能の意味

5. データの読み込みと期間の指定

import numpy as np
import xarray as xr
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
import cartopy.crs as ccrs
import cartopy.feature as cfeature

ncfile = "sst.mnmean.nc"
ds = xr.open_dataset(______)

sst = ds["sst"].values.astype(float)
lat = ds["lat"].values
lon = ds["lon"].values
time = pd.to_datetime(ds["time"].values)

sst[sst < -100] = np.nan

start_year = ______
end_year = ______

mask_period = (time.year >= start_year) & (time.year <= end_year)
sst_sub = sst[mask_period, :, :]
time_sub = time[mask_period]

ポイント

• xr.open_dataset で NetCDF を読む
• sst[sst < -100] = np.nan で欠損値を NaN にする
• 期間は start_year と end_year を変えるだけで指定できる

6. 年の小数への変換

回帰では時間を数値として扱いたいので、年月日ではなく「年の小数」に変換する。

years = time_sub.year + (time_sub.month - ______) / 12.0
years = years.to_numpy(dtype=float)

print(years[:5])
print(years[-5:])

ポイント

• 各月の中央付近を表すために 0.5 を使う
• たとえば 1993年1月は 1993.0416... のような値になる

7. 線形回帰と p 値の計算関数

ここでは、各格子点について x = year、y = SST とし、直線 y = ax + b をあてはめる。傾き a が SST トレンドであり、単位は °C/year である。

左：線形回帰によるトレンド計算の概念図。右：p 値（有意性）の概念図。

まず、何を検定するのか

この回の検定対象は、相関係数ではなく、回帰直線の傾きである。帰無仮説は次のように置く。

a = 0 は「時間が進んでも SST が増減しない」、つまりトレンドがないことを意味する。求めた傾きが 0 から十分離れているかを t 検定で調べる。

最小二乗法で傾き â を求める

観測データを (x_i, y_i) とする。ここでは x_i が年、y_i がその格子点の SST である。回帰直線による予測値を

と書く。観測値と予測値の差を残差という。

最小二乗法では、残差の二乗和

が最小になるように â と b̂ を決める。

傾き â の式が出てくるまで

まず、S を b̂ で微分し、最小値の条件として 0 とおく。

したがって、

両辺を n で割ると、

なので、切片は

となる。つまり、最小二乗法で求めた回帰直線は、平均点 (x̄, ȳ) を通る。

次に、S を â で微分して 0 とおく。

したがって、

ここに b̂ = ȳ − â x̄ を代入する。

括弧の中を整理すると、

ここで、

と置くと、微分条件は

である。一方、b̂ で微分した条件から、残差の和は 0 なので、

も成り立つ。したがって、

ここが、x_i を x_i − x̄ に置き換えてよい理由である。勝手に置き換えているのではなく、ΣA_i = 0 が成り立つため、同じ条件として書ける。

したがって、

これを展開すると、

よって、

傾きの標準誤差 SE(â)

傾き â を求めたら、次にその不確かさを見積もる。まず回帰直線から予測値を作る。

残差の二乗和は、

である。傾きと切片の 2 つをデータから推定しているため、自由度は

となる。残差分散は、

であり、傾きの標準誤差は、

である。点が直線のまわりに大きく散らばっているほど SSres が大きくなり、SE(â) も大きくなる。

傾きの t 検定

最後に、傾きが 0 からどれだけ離れているかを、傾きの不確かさで割って t 値を作る。

この式は、傾き â を求める式ではない。傾きはすでに最小二乗法で求めてあり、その傾きが 0 からどれだけはっきり離れているかを調べる式である。

Pythonコードとの対応

def calc_trend_stats(x, y):
    x = np.asarray(x, dtype=float)
    y = np.asarray(y, dtype=float)

    mask = np.isfinite(x) & np.isfinite(y)
    x = x[mask]
    y = y[mask]

    n = len(x)
    if n < 3:
        return np.nan, np.nan, np.nan, np.nan, np.nan

    # 最小二乗法で y = a*x + b を求める
    a, b = np.polyfit(x, y, 1)

    # 回帰直線による予測値
    y_pred = a * x + b

    # 残差平方和と全平方和
    ss_res = np.sum((y - y_pred) ** 2)
    ss_tot = np.sum((y - np.mean(y)) ** 2)
    r2 = np.nan if ss_tot == 0 else 1 - ss_res / ss_tot

    # x と y の相関係数
    r = np.corrcoef(x, y)[0, 1]

    # 傾き a の標準誤差と t 検定
    dof = n - 2
    sxx = np.sum((x - np.mean(x)) ** 2)
    mse = ss_res / dof
    se_a = np.sqrt(mse / sxx)
    tval = a / se_a
    pval = 2 * (1 - stats.t.cdf(np.abs(tval), df=dof))

    return r, a, b, r2, pval

ポイント

• a：傾き。単位は °C/year
• b：切片
• r2：決定係数
• pval：傾きが 0 かどうかを検定したときの p 値

8. 各格子点でのトレンド計算

nlat = len(lat)
nlon = len(lon)

trend = np.full((nlat, nlon), np.nan)
pval = np.full((nlat, nlon), np.nan)
rmap = np.full((nlat, nlon), np.nan)
r2map = np.full((nlat, nlon), np.nan)

min_count = ______

for i in range(nlat):
    for j in range(nlon):
        y = sst_sub[:, i, j]
        x = years

        mask = np.isfinite(y)
        if np.sum(mask) < min_count:
            continue

        r, a, b, r2, p = calc_trend_stats(x[mask], y[mask])

        trend[i, j] = ______
        pval[i, j] = ______
        rmap[i, j] = ______
        r2map[i, j] = ______

alpha = 0.05
sig_mask = pval < alpha

ポイント

• 外側のループは緯度、内側のループは経度
• 各格子で 1 本ずつ回帰直線を求める
• min_count は、有効データが少なすぎる格子を除くための条件である

9. トレンドマップの描画

fig = plt.figure(figsize=(14, 6))
ax = fig.add_subplot(1, 1, 1,
                     projection=ccrs.PlateCarree(central_longitude=______))

pcm = ax.pcolormesh(
    lon, lat, trend,
    transform=ccrs.PlateCarree(),
    shading="auto",
    cmap="coolwarm",
    vmin=-0.10, vmax=0.10
)

ax.add_feature(cfeature.LAND, facecolor="0.85", edgecolor="none", zorder=3)
ax.coastlines(linewidth=0.8, zorder=4)

gl = ax.gridlines(draw_labels=True, linewidth=0.5, color="0.7", alpha=0.8)
gl.top_labels = False
gl.right_labels = False

cbar = fig.colorbar(pcm, ax=ax, pad=0.02, shrink=0.9)
cbar.set_label("Trend (°C/year)")

ax.set_title(f"Global SST Trend (ERSSTv5, {start_year}-{end_year})")

fig.tight_layout()
fig.canvas.draw()
fig.savefig("sst_trend.png", dpi=300, bbox_inches="tight", facecolor="white")
plt.show()

ポイント

• central_longitude=180 にすると太平洋が中央にくる
• coolwarm では、暖色が正のトレンド、寒色が負のトレンドを表す
• Cartopy を使うときは、保存前に fig.canvas.draw() を入れると安全である

10. p 値と有意マスクの表示

有意な格子点だけを黒点で重ねて表示する。

yy, xx = np.where(sig_mask)

ax.scatter(
    lon[xx], lat[yy],
    s=2, c="k", marker="o",
    transform=ccrs.PlateCarree(),
    zorder=5, alpha=_____
)

ポイント

• np.where(sig_mask) は True の場所の添字を返す
• 黒点は「p 値が十分小さい格子点」を表す
• ただし、p 値が小さいことと、物理的に重要であることは同じではない

11. p 値とは何か

今回の検定では、各格子点について次の帰無仮説を考えている。

つまり、「その格子点では長期トレンドがない」と仮定する。この仮定のもとで、実際に得られた傾きと同じくらい、またはそれ以上に 0 から離れた傾きが偶然に出る確率が p 値である。

• p 値が小さい：傾き 0 では説明しにくい。トレンドがあると考えやすい。
• p 値が大きい：傾き 0 でも十分あり得る。トレンドがあるとは言いにくい。

たとえば alpha = 0.05 としたとき、p < 0.05 なら「5%水準で有意」と呼ぶことが多い。

12. 議論：SST変動は何で説明できるか

ここまでで、各格子点における海面水温（SST）の長期トレンドを求め、 全球平均が上昇していても、その増減の大きさや有意性は空間的に一様ではないことを確認した。 では、この空間分布の違いは何によって説明できるのだろうか。

短期変動としての ENSO

前回の ENSO コンポジット解析では、El Niño と La Niña のときに 海面水温の分布がどのように変わるかを調べた。 これは、SST の年々変動、つまり比較的短い時間スケールの変動を見ている。

たとえば赤道太平洋では ENSO の影響が強く、ある年には暖かく、 別の年には冷たくなる。このような場所では、SST の分布を理解するうえで ENSO が重要な説明要因になる。

長期変化としてのトレンド

一方、このページで見ているトレンドは、数十年スケールで SST がどの方向に変化してきたかを表している。 これは、温暖化や海洋循環の長期変化など、 長い時間スケールの変動を反映している可能性がある。

したがって、トレンドマップは「この海域では長い目で見ると暖まりやすいのか、 それともあまり変化していないのか」を示している。

ENSOだけで説明できるか

ENSO コンポジットで大きな差が見られる海域では、 年々変動のかなりの部分を ENSO が説明していると考えられる。 しかし、ENSO の影響が弱い海域でも長期トレンドが強く現れることがある。 その場合、SST の変化は ENSO だけでは説明できない。

トレンドだけで説明できるか

逆に、長期トレンドが見られても、年ごとのばらつきが非常に大きい海域では、 実際の SST 変動をトレンドだけで説明することは難しい。 たとえば ENSO のような大きな年々変動が重なると、 長期的には上昇傾向でも、ある年には一時的に低温になることがある。

この回で得たい理解

• ENSO は SST の年々変動を説明する重要な要因である
• トレンドは SST の長期変化を説明する重要な要因である
• 海域によって、どちらが効いているかは異なる
• 一つの要因だけでは、全球の SST 分布を説明できない

このように、SST の変動を理解するには、 短期変動（ENSO）と長期変化（トレンド）を分けて考えることが重要である。 今後はさらに、海面高度や流れの情報も使いながら、 「なぜその場所で変化が起きるのか」を物理的に考えていく。

13. まとめ

• ERSSTv5 の月平均データを NetCDF から読み込んだ
• 時間を年の小数に変換し、各格子点で時間方向の線形回帰を行った
• 傾き a を SST トレンドとして地図に表した
• p 値を計算し、統計的に有意な格子点を重ね描きした
• 全球平均が上昇していても、空間分布は一様ではないことを確認した

14. 解答の表示

まずは上の穴埋めコードを自分で考えること。確認したい場合だけ、パスワードを入力して解答と完全スクリプトを表示する。

ncfile
1993
2021
0.5
30
a
p
r
r2
180
0.6

import numpy as np
import xarray as xr
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
import cartopy.crs as ccrs
import cartopy.feature as cfeature

ncfile = "sst.mnmean.nc"
ds = xr.open_dataset(ncfile)

sst = ds["sst"].values.astype(float)
lat = ds["lat"].values
lon = ds["lon"].values
time = pd.to_datetime(ds["time"].values)

sst[sst < -100] = np.nan

start_year = 1993
end_year = 2021

mask_period = (time.year >= start_year) & (time.year <= end_year)
sst_sub = sst[mask_period, :, :]
time_sub = time[mask_period]

years = time_sub.year + (time_sub.month - 0.5) / 12.0
years = years.to_numpy(dtype=float)

def calc_trend_stats(x, y):
    x = np.asarray(x, dtype=float)
    y = np.asarray(y, dtype=float)

    mask = np.isfinite(x) & np.isfinite(y)
    x = x[mask]
    y = y[mask]

    n = len(x)
    if n < 3:
        return np.nan, np.nan, np.nan, np.nan, np.nan

    a, b = np.polyfit(x, y, 1)

    y_pred = a * x + b
    ss_res = np.sum((y - y_pred) ** 2)
    ss_tot = np.sum((y - np.mean(y)) ** 2)
    r2 = np.nan if ss_tot == 0 else 1 - ss_res / ss_tot

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

    dof = n - 2
    sxx = np.sum((x - np.mean(x)) ** 2)

    mse = ss_res / dof
    se_a = np.sqrt(mse / sxx)
    tval = a / se_a
    pval = 2 * (1 - stats.t.cdf(np.abs(tval), df=dof))

    return r, a, b, r2, pval

nlat = len(lat)
nlon = len(lon)

trend = np.full((nlat, nlon), np.nan)
pval = np.full((nlat, nlon), np.nan)
rmap = np.full((nlat, nlon), np.nan)
r2map = np.full((nlat, nlon), np.nan)

min_count = 30

for i in range(nlat):
    for j in range(nlon):
        y = sst_sub[:, i, j]
        x = years

        mask = np.isfinite(y)
        if np.sum(mask) < min_count:
            continue

        r, a, b, r2, p = calc_trend_stats(x[mask], y[mask])

        trend[i, j] = a
        pval[i, j] = p
        rmap[i, j] = r
        r2map[i, j] = r2

alpha = 0.05
sig_mask = pval < alpha

fig = plt.figure(figsize=(14, 6))
ax = fig.add_subplot(1, 1, 1,
                     projection=ccrs.PlateCarree(central_longitude=180))

pcm = ax.pcolormesh(
    lon, lat, trend,
    transform=ccrs.PlateCarree(),
    shading="auto",
    cmap="coolwarm",
    vmin=-0.10, vmax=0.10
)

ax.add_feature(cfeature.LAND, facecolor="0.85", edgecolor="none", zorder=3)
ax.coastlines(linewidth=0.8, zorder=4)

gl = ax.gridlines(draw_labels=True, linewidth=0.5, color="0.7", alpha=0.8)
gl.top_labels = False
gl.right_labels = False

cbar = fig.colorbar(pcm, ax=ax, pad=0.02, shrink=0.9)
cbar.set_label("Trend (°C/year)")

yy, xx = np.where(sig_mask)
ax.scatter(
    lon[xx], lat[yy],
    s=2, c="k", marker="o",
    transform=ccrs.PlateCarree(),
    zorder=5, alpha=0.6
)

ax.set_title(f"Global SST Trend (ERSSTv5, {start_year}-{end_year})\nblack dots: p < {alpha}")

fig.tight_layout()
fig.canvas.draw()
fig.savefig("sst_trend.png", dpi=300, bbox_inches="tight", facecolor="white")
plt.show()