データサイエンス 2026 / 11 SLAから地衡流・相対渦度・EKEを計算する

11 SLAから地衡流・相対渦度・EKEを計算する

海面高度偏差（SLA）から地衡流を求め、黒潮・黒潮続流・東北沖の暖水塊にともなう相対渦度と渦運動エネルギー（EKE）を可視化する。

講義名
データサイエンス
使用言語
Python / xarray / cartopy
到達目標
SLAから流速・渦度・EKEを計算し、海洋循環の力学的特徴を図から説明できるようになる

1. 背景：SLAから何がわかるか

前回までは、海面高度偏差（SLA）を使ってENSOや長期トレンドを調べた。 今回はSLAの空間勾配から地衡流を計算し、海の流れの構造を可視化する。

海面高度が高い場所と低い場所があると、圧力傾度力が生じる。地球が回転しているため、 大規模な海洋ではこの圧力傾度力とコリオリ力がほぼ釣り合い、地衡流が生じる。

今回の狙いは、黒潮、黒潮続流、東北沖の暖水塊を「高さの偏差」だけでなく、「流れ」「回転」「渦活動」として見ることである。

2. 地衡流・相対渦度・EKEの式

SLAを η、重力加速度を g、コリオリパラメータを f とすると、SLA由来の地衡流は次のように書ける。

ug = −(g/f) ∂η/∂y
vg = (g/f) ∂η/∂x

ここで、ug は東向き流速、vg は北向き流速である。

ζ = ∂vg/∂x − ∂ug/∂y

ζ は相対渦度で、流れの回転の強さを表す。

EKE = 1/2 (ug2 + vg2)
ここではSLAから計算するため、厳密には「SLA由来の偏差地衡流」と「その運動エネルギー」である。平均流まで含めたい場合は、SLAではなくADTを用いる。

3. 使用データと作成する図

このページと同じディレクトリに次のファイルを置く。

SLAデータ (sla_monthly_0125deg_2024.nc)
ファイル内容
sla_monthly_0125deg_2024.nc2024年1月〜12月の月平均SLA。空間解像度は0.125度。

この課題では以下の4枚を作成する。

  1. 2024年5月のSLAとSLA由来地衡流
  2. 2024年5月の相対渦度
  3. 2024年5月のEKE
  4. 2024年平均EKE
STEP 1

4. ライブラリと設定

import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature

infile = "______________________________.nc"
outfile = "sla_geostrophy_vorticity_eke_japan_2024.nc"

lon_min, lon_max = ______, ______
lat_min, lat_max = ______, ______

month_to_plot = ______
g = ______
omega = ______
Re = ______

ポイント

STEP 2

5. SLAデータの読み込み

ds = xr.open_dataset(infile)
sla = ds["____"].astype("float64")

sla = sla.sel(
    longitude=slice(lon_min-2, lon_max+2),
    latitude=slice(lat_min-2, lat_max+2)
)

lon = sla["_________"]
lat = sla["________"]

print(sla)

ポイント

STEP 3

6. 緯度経度を距離に換算する

lat_rad = np.deg2rad(lat)

f = 2 * omega * np.sin(________)
f = xr.DataArray(f, coords={"latitude": lat}, dims=("latitude",))

meters_per_deg_lat = np.pi * Re / 180.0
meters_per_deg_lon = meters_per_deg_lat * np.cos(________)
meters_per_deg_lon = xr.DataArray(
    meters_per_deg_lon,
    coords={"latitude": lat},
    dims=("latitude",)
)

ポイント

STEP 4

7. SLA勾配から地衡流を計算する

deta_dlat = sla.differentiate("________") / meters_per_deg_lat
deta_dlon = sla.differentiate("_________") / meters_per_deg_lon

ug = -(g / f) * __________
vg =  (g / f) * __________

ug.name = "ug"
vg.name = "vg"
ug.attrs["units"] = "m s-1"
vg.attrs["units"] = "m s-1"

ポイント

STEP 5

8. 相対渦度とEKEを計算する

dvdx = vg.differentiate("longitude") / meters_per_deg_lon
dudy = ug.differentiate("latitude") / meters_per_deg_lat

zeta = ______ - ______
zeta.name = "relative_vorticity"
zeta.attrs["units"] = "s-1"

zeta_over_f = zeta / ____
zeta_over_f.name = "relative_vorticity_over_f"

eke = 0.5 * ((ug * ______)**2 + (vg * ______)**2)
eke.name = "eke"
eke.attrs["units"] = "cm2 s-2"

log10_mean_eke = np.log10(eke.mean("____"))

ポイント

STEP 6

9. NetCDFとして保存する

sla_p = sla.sel(longitude=slice(lon_min, lon_max), latitude=slice(lat_min, lat_max))
ug_p = ug.sel(longitude=slice(lon_min, lon_max), latitude=slice(lat_min, lat_max))
vg_p = vg.sel(longitude=slice(lon_min, lon_max), latitude=slice(lat_min, lat_max))
zeta_p = zeta.sel(longitude=slice(lon_min, lon_max), latitude=slice(lat_min, lat_max))
eke_p = eke.sel(longitude=slice(lon_min, lon_max), latitude=slice(lat_min, lat_max))
log10_mean_eke_p = log10_mean_eke.sel(longitude=slice(lon_min, lon_max), latitude=slice(lat_min, lat_max))

out = xr.Dataset({
    "sla": sla_p.astype("float32"),
    "ug": ug_p.astype("float32"),
    "vg": vg_p.astype("float32"),
    "relative_vorticity": zeta_p.astype("float32"),
    "relative_vorticity_over_f": zeta_over_f.sel(
        longitude=slice(lon_min, lon_max),
        latitude=slice(lat_min, lat_max)
    ).astype("float32"),
    "eke": eke_p.astype("float32"),
    "log10_mean_eke": log10_mean_eke_p.astype("float32"),
})

out.to_netcdf(________)
print(f"Saved: {outfile}")
STEP 7

10. 日本周辺で描画する

def plot_map(field, u, v, title, cbar_label,
             vmin=None, vmax=None,
             cmap="turbo", outpng="figure.png",
             vector_step=8,
             contour_sla=None):

    proj = ccrs.PlateCarree()
    fig = plt.figure(figsize=(10, 8))
    ax = plt.axes(projection=proj)
    ax.set_extent([lon_min, lon_max, lat_min, lat_max], crs=proj)

    im = ax.pcolormesh(
        field["longitude"],
        field["latitude"],
        field,
        transform=proj,
        cmap=cmap,
        vmin=vmin,
        vmax=vmax,
        shading="auto"
    )

    ax.coastlines(resolution="10m", linewidth=0.8, zorder=20)
    ax.add_feature(cfeature.LAND, facecolor="lightgray", edgecolor="black", linewidth=0.5, zorder=10)

    if contour_sla is not None:
        cs = ax.contour(
            contour_sla["longitude"],
            contour_sla["latitude"],
            contour_sla,
            levels=np.arange(-1.0, 1.05, 0.1),
            colors="k",
            linewidths=0.45,
            alpha=0.55,
            transform=proj
        )

    uu = u.isel(latitude=slice(None, None, vector_step),
                longitude=slice(None, None, vector_step))
    vv = v.isel(latitude=slice(None, None, vector_step),
                longitude=slice(None, None, vector_step))

    speed = np.sqrt(uu**2 + vv**2)
    mask = speed > ______

    q = ax.quiver(
        uu["longitude"], uu["latitude"],
        uu.where(mask), vv.where(mask),
        transform=proj,
        color="k",
        scale=13,
        width=0.0014,
        alpha=0.70,
        zorder=30
    )

    ax.quiverkey(
        q,
        X=0.08, Y=0.87,
        U=______,
        label="0.5 m/s",
        labelpos="E",
        coordinates="axes",
        fontproperties={"size": 12, "weight": "bold"}
    )

    gl = ax.gridlines(draw_labels=True, linewidth=0.5, alpha=0.45)
    gl.top_labels = False
    gl.right_labels = False

    ax.set_title(title, fontsize=18, weight="bold")
    cb = plt.colorbar(im, ax=ax, orientation="vertical", pad=0.02)
    cb.set_label(cbar_label)

    plt.tight_layout()
    plt.savefig(outpng, dpi=300, bbox_inches="tight", facecolor="white")
    plt.show()

ポイント

作成例

SLA-derived geostrophic velocity
2024年5月のSLAとSLA由来地衡流。黒潮続流と暖水塊にともなう流れが見える。
Relative vorticity
2024年5月の相対渦度。暖水塊・冷水塊の周辺で正負の渦度が強く現れる。
EKE
2024年5月のEKE。黒潮・黒潮続流域で渦運動エネルギーが大きい。
Mean EKE
2024年平均EKE。月ごとのノイズがならされ、黒潮続流域の高EKE帯が明瞭になる。

11. 考察課題

考察せよ
  • 黒潮・黒潮続流は、SLAの高低差とどのように対応しているか。
  • 東北沖の暖水塊では、SLA、流速ベクトル、相対渦度はどのような関係になっているか。
  • 相対渦度図で、正の渦度と負の渦度はどのような場所に現れているか。
  • 月平均EKEと年平均EKEでは、どちらが黒潮続流域の渦活動帯を説明しやすいか。
  • SLA由来の地衡流と、ADT由来の地衡流では何が違うと考えられるか。
SLA由来の地衡流は、平均流そのものではなく偏差流を表す。黒潮の平均的な流路をより直接的に見たい場合は、ADTを用いた同じ計算を行う必要がある。

12. 解答の表示

課題提出後、確認用としてパスワードを入力すると、穴埋めの要点と完全スクリプトを表示する。