11 黒潮海域の海面高度と地衡流 | データサイエンス 2026

1. 目的：実習前に黒潮を読む

海洋学実習では、船で観測する前に、衛星データから現在の海況を把握しておくことが重要である。この回では、黒潮海域の海面高度データを使い、黒潮の位置、海面高度勾配、地衡流を自分で確認する。

このページの到達目標は、海面高度から地衡流を計算し、CMEMSに含まれる流速と比較できるようになることである。

2. 使用データ

使用するファイルは、作業ディレクトリ直下に置く。

kuroshio_altimetry_daily_20260501_20260512.nc 参考PDF：地衡流計算のイメージ

変数	意味	今回の使い方
`adt`	Absolute Dynamic Topography。平均場を含む海面高度。	黒潮の流路と絶対地衡流を見る。
`sla`	Sea Level Anomaly。平均海面からの偏差。	渦・蛇行・地衡流偏差を見る。
`ugos, vgos`	CMEMSに含まれる絶対地衡流。	ADTから自分で計算した流速と比較する。
`ugosa, vgosa`	CMEMSに含まれる地衡流偏差。	SLAから自分で計算した流速と比較する。

ファイルの置き場所：この回では、学生がパスで混乱しないように、NetCDFファイルをサブディレクトリではなく、ノートブックと同じ作業ディレクトリに置く。

3. 海面高度から地衡流へ

海面高度が場所によって違うと、海中には水平圧力勾配が生じる。静水圧近似では、鉛直方向の圧力勾配は重力とつり合うため、海面の傾きは水平方向の圧力勾配として現れる。

静水圧近似： ∂p/∂z ≃ −ρg

海面の傾きによる水平圧力勾配：
∂p/∂x ≃ ρg ∂η/∂x, ∂p/∂y ≃ ρg ∂η/∂y

時間変化項や非線形項がコリオリ項に比べて小さい場合、水平運動は主に圧力勾配力とコリオリ力のつり合いで近似できる。この近似が地衡流である。

u_g = − (g/f) ∂η/∂y, v_g = (g/f) ∂η/∂x

ここで、u_g は東西方向流速、v_g は南北方向流速、g は重力加速度、f = 2Ωsinφ はコリオリパラメータ、η は海面高度である。

北半球では、流体はおおまかに「高い海面を右手に見て」流れる。黒潮は、海面高度の高い亜熱帯側と低い亜寒帯側の境界付近に現れる強い流れとして見える。

4. 有限差分：前方差分と中央差分

コンピュータ上の海面高度は、連続関数ではなく格子点上の値である。そのため、偏微分をそのまま計算するのではなく、格子点間の差で近似する。これを有限差分という。

Forward and central difference in one dimension

図1. 1次元で見た前方差分と中央差分。前方差分は注目点と隣の点、中央差分は注目点の両側の点を使う。

Grid points used for one-sided forward difference

図2. 2次元格子での前方差分。注目点と右隣・上隣を使って勾配を近似する。

図3. 2次元格子での中央差分。左右・上下の格子点を使うため、外側1格子は計算できない。

用語の注意：ここで比較する onesided は、厳密な意味での移流方程式の風上差分ではなく、海面高度勾配を片側の格子差で近似する片側差分である。授業では、中央差分との違いを見るための比較用として使う。

5. 実行例：今回使う図

ページ内では、まずADTとSLAの基本図を見せ、そのあとCMEMS流速と自分で計算した地衡流の重ね描き、さらに差分図を示す。本文の中心に使うのは図4〜図8で、図9〜図10は差分法の違いを考える発展例として扱う。

図4. ADT（Absolute Dynamic Topography）。黒潮は海面高度の等値線が密な帯として見える。

図5. SLA（Sea Level Anomaly）。平均場からの偏差として、黒潮蛇行や渦構造が見える。

図6. ADTに絶対地衡流を重ねた図。黒矢印はCMEMS流速、赤矢印はADTから計算した流速。

図7. SLAに地衡流偏差を重ねた図。SLAから計算した偏差流速とCMEMSの偏差流速を比較する。

図8. ADTから計算した東西流速とCMEMS ugos の差。黒潮流軸や沿岸で差が目立つ。

図9. ADTで片側差分と中央差分を比較した例。強い勾配域では差分法の違いが流速差として現れる。

図10. SLAで片側差分と中央差分を比較した例。小さな差でも、格子スケールの数値微分の違いが見える。

6. 完全スクリプト

下のコードをJupyter Labのセルにコピーして実行する。図は画面に表示されるだけでなく、figures_11_kuroshio/ フォルダにもPNGとして保存される。

注意：このスクリプトは、kuroshio_altimetry_daily_20260501_20260512.nc がノートブックと同じ作業ディレクトリにあることを前提にしている。

# ============================================================
# 11. 黒潮海域の海面高度・地衡流解析（完全スクリプト）
# ============================================================
#
# このスクリプトで行うこと
# ------------------------------------------------------------
# 1. 黒潮海域の海面高度データ（NetCDF）を読み込む
# 2. ADT（Absolute Dynamic Topography）を描く
# 3. SLA（Sea Level Anomaly）を描く
# 4. CMEMSの流速（ugos/vgos, ugosa/vgosa）を表示する
# 5. ADT / SLA から自分で地衡流を計算する
# 6. 差分法を中央差分 / 片側差分で切り替える
# 7. CMEMS流速と自分で計算した流速を比較する
# 8. 差の統計量や差分マップを描く
# 9. Webサイト掲載用に各図をPNG保存する
#
# 注意
# ------------------------------------------------------------
# ・厳密には "風上差分" ではなく、ここでは "片側差分" を比較用として使う。
# ・教育上、「差分法によって数値微分の結果が変わる」ことを見せるのが目的。
# ・レジェンドは左上（中国大陸側）に配置する。
# ============================================================


# ============================================================
# 0. ライブラリの読み込み
# ============================================================

import numpy as np
import xarray as xr
import matplotlib.pyplot as plt

import cartopy.crs as ccrs
import cartopy.feature as cfeature

from matplotlib.colors import LinearSegmentedColormap, TwoSlopeNorm
from matplotlib.lines import Line2D
from pathlib import Path


# ============================================================
# 1. 基本設定
# ============================================================

# ファイル名（作業ディレクトリ直下に置く）
filename = "kuroshio_altimetry_daily_20260501_20260512.nc"

# 使う日
# it = 0 で最初の日、it = -1 で最後の日
it = -1

# 地図範囲
lon_min, lon_max = 100, 150
lat_min, lat_max = 20, 50

# ベクトルの間引き
skip = 16

# ベクトルの長さ調整
# scale を大きくすると矢印は短くなる
vector_scale = 18

# 比較に使う差分法
# "central"   : 中央差分
# "onesided"  : 片側差分（比較用）
method_main = "central"

# 片側差分と中央差分の違いも見たい場合
compare_two_methods = True

# 図の保存先フォルダ
figdir = Path("figures_11_kuroshio")
figdir.mkdir(exist_ok=True)

# 図保存用の通し番号
fig_count = 0

def save_current_figure(label):
    """
    現在の図をPNGとして保存する。
    plt.show() の直前で呼び出す。
    """
    global fig_count

    fig_count += 1
    out_png = figdir / f"fig11_kuroshio_{fig_count:02d}_{label}_{date_str}.png"

    plt.gcf().savefig(
        out_png,
        dpi=300,
        bbox_inches="tight",
        facecolor="white"
    )

    print("Saved figure:", out_png)


# 定数
g = 9.81             # 重力加速度 [m/s^2]
Omega = 7.2921e-5    # 地球自転角速度 [rad/s]
R = 6371000.0        # 地球半径 [m]


# ============================================================
# 2. 海面高度らしいカラーマップ
# ============================================================
#
# ADTは「海面高度らしい」色：
#   低い = 青
#   高い = 黄～赤
#
# SLAは偏差なので、0を中心に
#   負 = 青
#   正 = 赤
# ============================================================

ssh_colors = [
    "#08306b",  # deep blue
    "#2171b5",  # blue
    "#6baed6",  # light blue
    "#c7e9f1",  # very light cyan
    "#ffffbf",  # pale yellow
    "#fdae61",  # orange
    "#d7191c"   # red
]
ssh_cmap = LinearSegmentedColormap.from_list("ssh_cmap", ssh_colors, N=256)
sla_cmap = "RdBu_r"


# ============================================================
# 3. データ読み込み
# ============================================================

ds = xr.open_dataset(filename)

print(ds)

lon = ds["longitude"].values
lat = ds["latitude"].values

date_str = str(ds["time"].values[it])[:10]
print("Selected date:", date_str)

adt   = ds["adt"].isel(time=it)
sla   = ds["sla"].isel(time=it)
ugos  = ds["ugos"].isel(time=it)
vgos  = ds["vgos"].isel(time=it)
ugosa = ds["ugosa"].isel(time=it)
vgosa = ds["vgosa"].isel(time=it)


# ============================================================
# 4. 地図の共通設定関数
# ============================================================

def setup_map(ax):
    """
    地図の共通設定
    """
    ax.set_extent([lon_min, lon_max, lat_min, lat_max], crs=ccrs.PlateCarree())

    # 陸を灰色で塗る
    ax.add_feature(cfeature.LAND, facecolor="0.82", zorder=5)

    # 海岸線
    ax.coastlines(resolution="50m", linewidth=0.9, zorder=6)

    # 緯度経度線
    gl = ax.gridlines(
        draw_labels=True,
        linewidth=0.5,
        color="gray",
        alpha=0.5,
        linestyle="--"
    )
    gl.top_labels = False
    gl.right_labels = False

    return ax


# ============================================================
# 5. 差分法を切り替えられる地衡流計算関数
# ============================================================

def calc_geostrophic_current_switch(eta, lon, lat, method="central"):
    eta = np.asarray(eta, dtype=float)

    # 勾配配列を NaN で初期化
    d_eta_dx = np.full_like(eta, np.nan, dtype=float)
    d_eta_dy = np.full_like(eta, np.nan, dtype=float)

    # 緯度経度をラジアンに変換
    lon_rad = np.deg2rad(lon)
    lat_rad = np.deg2rad(lat)

    lon2d, lat2d = np.meshgrid(lon, lat)
    lat2d_rad = np.deg2rad(lat2d)

    # コリオリパラメータ
    f = 2.0 * Omega * np.sin(lat2d_rad)

    # 有効点判定
    valid = np.isfinite(eta)

    # --------------------------------------------------------
    # 中央差分
    # --------------------------------------------------------
    if method == "central":

        # x方向勾配 dη/dx
        # [η(j+1) - η(j-1)] / [x(j+1)-x(j-1)]
        dlon2 = lon_rad[2:] - lon_rad[:-2]
        dx2d = R * np.cos(lat2d_rad[:, 1:-1]) * dlon2[np.newaxis, :]
        d_eta_dx[:, 1:-1] = (eta[:, 2:] - eta[:, :-2]) / dx2d

        # y方向勾配 dη/dy
        # [η(i+1) - η(i-1)] / [y(i+1)-y(i-1)]
        dlat2 = lat_rad[2:] - lat_rad[:-2]
        dy2d = R * dlat2[:, np.newaxis]
        d_eta_dy[1:-1, :] = (eta[2:, :] - eta[:-2, :]) / dy2d

        # 周囲4方向の値が必要
        valid_grad = np.zeros_like(valid, dtype=bool)
        valid_grad[1:-1, 1:-1] = (
            valid[1:-1, 1:-1] &
            valid[1:-1, :-2] &
            valid[1:-1, 2:] &
            valid[:-2, 1:-1] &
            valid[2:, 1:-1]
        )

    # --------------------------------------------------------
    # 片側差分（比較用）
    # --------------------------------------------------------
    elif method == "onesided":

        # x方向勾配 dη/dx
        # [η(j+1) - η(j)] / [x(j+1)-x(j)]
        dlon1 = lon_rad[1:] - lon_rad[:-1]
        dx2d = R * np.cos(lat2d_rad[:, :-1]) * dlon1[np.newaxis, :]
        d_eta_dx[:, :-1] = (eta[:, 1:] - eta[:, :-1]) / dx2d

        # y方向勾配 dη/dy
        # [η(i+1) - η(i)] / [y(i+1)-y(i)]
        dlat1 = lat_rad[1:] - lat_rad[:-1]
        dy2d = R * dlat1[:, np.newaxis]
        d_eta_dy[:-1, :] = (eta[1:, :] - eta[:-1, :]) / dy2d

        # 右隣・上隣の値が必要
        valid_grad = np.zeros_like(valid, dtype=bool)
        valid_grad[:-1, :-1] = (
            valid[:-1, :-1] &
            valid[:-1, 1:] &
            valid[1:, :-1]
        )

    else:
        raise ValueError('method must be "central" or "onesided"')

    # 地衡流
    ug = -g / f * d_eta_dy
    vg =  g / f * d_eta_dx

    # 無効点はNaN
    ug[~valid_grad] = np.nan
    vg[~valid_grad] = np.nan

    speed = np.sqrt(ug**2 + vg**2)

    return ug, vg, speed


# ============================================================
# 6. 比較統計を出す関数
# ============================================================

def compare_velocity(u_calc, v_calc, u_ref, v_ref, name):
    """
    計算した流速と参照流速（CMEMS）を比較する。
    """
    valid = (
        np.isfinite(u_calc) &
        np.isfinite(v_calc) &
        np.isfinite(u_ref) &
        np.isfinite(v_ref)
    )

    du = np.full_like(u_calc, np.nan)
    dv = np.full_like(v_calc, np.nan)

    du[valid] = u_calc[valid] - u_ref[valid]
    dv[valid] = v_calc[valid] - v_ref[valid]

    speed_calc = np.sqrt(u_calc**2 + v_calc**2)
    speed_ref  = np.sqrt(u_ref**2 + v_ref**2)

    dspeed = np.full_like(speed_calc, np.nan)
    dspeed[valid] = speed_calc[valid] - speed_ref[valid]

    print("============================================================")
    print(name)
    print("============================================================")
    print("Number of valid points:", np.sum(valid))
    print("du mean   :", np.nanmean(du[valid]), "m/s")
    print("dv mean   :", np.nanmean(dv[valid]), "m/s")
    print("du RMSE   :", np.sqrt(np.nanmean(du[valid]**2)), "m/s")
    print("dv RMSE   :", np.sqrt(np.nanmean(dv[valid]**2)), "m/s")
    print("speed RMSE:", np.sqrt(np.nanmean(dspeed[valid]**2)), "m/s")
    print("du max abs:", np.nanmax(np.abs(du[valid])), "m/s")
    print("dv max abs:", np.nanmax(np.abs(dv[valid])), "m/s")

    corr_u = np.corrcoef(u_calc[valid], u_ref[valid])[0, 1]
    corr_v = np.corrcoef(v_calc[valid], v_ref[valid])[0, 1]

    print("u corr    :", corr_u)
    print("v corr    :", corr_v)

    return du, dv, dspeed


# ============================================================
# 7. 差分マップ描画関数
# ============================================================

def plot_difference_map(diff, title, label, vmax=None):
    """
    差分マップを描く。
    """
    fig = plt.figure(figsize=(10, 7))
    ax = plt.axes(projection=ccrs.PlateCarree())
    setup_map(ax)

    if vmax is None:
        vmax = float(np.nanpercentile(np.abs(diff), 98))

    p = ax.pcolormesh(
        lon, lat, diff,
        transform=ccrs.PlateCarree(),
        shading="auto",
        cmap="RdBu_r",
        norm=TwoSlopeNorm(vmin=-vmax, vcenter=0.0, vmax=vmax)
    )

    cb = plt.colorbar(p, ax=ax, shrink=0.85, pad=0.03)
    cb.set_label(label)

    plt.title(title)

    save_current_figure("difference")

    plt.show()


# ============================================================
# 8. ベクトル重ね描き関数
# ============================================================

def overlay_two_vector_fields(
    background,
    u1, v1,
    u2, v2,
    title,
    cbar_label,
    cmap,
    contour_levels=None,
    vmin=None,
    vmax=None,
    norm=None
):
    fig = plt.figure(figsize=(10, 7))
    ax = plt.axes(projection=ccrs.PlateCarree())
    setup_map(ax)

    # 背景図
    if norm is None:
        p = ax.pcolormesh(
            lon, lat, background,
            transform=ccrs.PlateCarree(),
            shading="auto",
            cmap=cmap,
            vmin=vmin,
            vmax=vmax
        )
    else:
        p = ax.pcolormesh(
            lon, lat, background,
            transform=ccrs.PlateCarree(),
            shading="auto",
            cmap=cmap,
            norm=norm
        )

    # 等値線
    if contour_levels is not None:
        cs = ax.contour(
            lon, lat, background,
            levels=contour_levels,
            colors="0.35",
            linewidths=0.6,
            transform=ccrs.PlateCarree(),
            zorder=4
        )
        ax.clabel(cs, fontsize=8, fmt="%.2f")

    # ベクトルを間引く
    lon_q = lon[::skip]
    lat_q = lat[::skip]

    u1_q = u1[::skip, ::skip]
    v1_q = v1[::skip, ::skip]

    u2_q = u2[::skip, ::skip]
    v2_q = v2[::skip, ::skip]

    # 黒：CMEMS
    q1 = ax.quiver(
        lon_q, lat_q,
        u1_q, v1_q,
        transform=ccrs.PlateCarree(),
        color="black",
        scale=vector_scale,
        width=0.0022,
        headwidth=3.8,
        headlength=4.8,
        zorder=7
    )

    # 赤：自分で計算した流速
    q2 = ax.quiver(
        lon_q, lat_q,
        u2_q, v2_q,
        transform=ccrs.PlateCarree(),
        color="red",
        scale=vector_scale,
        width=0.0020,
        headwidth=3.8,
        headlength=4.8,
        alpha=0.9,
        zorder=8
    )

    # 凡例（左上：中国大陸側）
    legend_elements = [
        Line2D([0], [0], color="black", lw=2, label="CMEMS velocity"),
        Line2D([0], [0], color="red",   lw=2, label="Calculated from height")
    ]

    ax.legend(
        handles=legend_elements,
        loc="upper left",
        bbox_to_anchor=(0.02, 0.98),
        frameon=True,
        framealpha=0.95,
        fontsize=11
    )

    # 矢印スケール
    ax.quiverkey(
        q1,
        X=0.88, Y=-0.08,
        U=1.0,
        label="1.0 m/s",
        labelpos="E",
        coordinates="axes"
    )

    cb = plt.colorbar(p, ax=ax, shrink=0.85, pad=0.03)
    cb.set_label(cbar_label)

    plt.title(title)

    save_current_figure("overlay")

    plt.show()


# ============================================================
# 9. ADT のマッピング
# ============================================================

adt_vmin = float(np.nanpercentile(adt, 2))
adt_vmax = float(np.nanpercentile(adt, 98))

fig = plt.figure(figsize=(10, 7))
ax = plt.axes(projection=ccrs.PlateCarree())
setup_map(ax)

p = ax.pcolormesh(
    lon, lat, adt,
    transform=ccrs.PlateCarree(),
    shading="auto",
    cmap=ssh_cmap,
    vmin=adt_vmin,
    vmax=adt_vmax
)

cs = ax.contour(
    lon, lat, adt,
    levels=np.arange(0.2, 2.1, 0.1),
    colors="black",
    linewidths=0.5,
    transform=ccrs.PlateCarree()
)
ax.clabel(cs, fontsize=8, fmt="%.1f")

cb = plt.colorbar(p, ax=ax, shrink=0.85, pad=0.03)
cb.set_label("ADT (m)")

plt.title(f"Absolute Dynamic Topography (ADT)\n{date_str}")

save_current_figure("adt")

plt.show()


# ============================================================
# 10. SLA のマッピング
# ============================================================

sla_abs = float(np.nanpercentile(np.abs(sla), 98))

fig = plt.figure(figsize=(10, 7))
ax = plt.axes(projection=ccrs.PlateCarree())
setup_map(ax)

p = ax.pcolormesh(
    lon, lat, sla,
    transform=ccrs.PlateCarree(),
    shading="auto",
    cmap=sla_cmap,
    norm=TwoSlopeNorm(vmin=-sla_abs, vcenter=0.0, vmax=sla_abs)
)

cs = ax.contour(
    lon, lat, sla,
    levels=np.arange(-0.5, 0.55, 0.05),
    colors="black",
    linewidths=0.4,
    transform=ccrs.PlateCarree()
)
ax.clabel(cs, fontsize=8, fmt="%.2f")

cb = plt.colorbar(p, ax=ax, shrink=0.85, pad=0.03)
cb.set_label("SLA (m)")

plt.title(f"Sea Level Anomaly (SLA)\n{date_str}")

save_current_figure("sla")

plt.show()


# ============================================================
# 11. ADT から地衡流を計算
# ============================================================

ug_adt_main, vg_adt_main, speed_adt_main = calc_geostrophic_current_switch(
    adt.values, lon, lat, method=method_main
)

speed_cmems_abs = np.sqrt(ugos.values**2 + vgos.values**2)


# ============================================================
# 12. SLA から地衡流偏差を計算
# ============================================================

ug_sla_main, vg_sla_main, speed_sla_main = calc_geostrophic_current_switch(
    sla.values, lon, lat, method=method_main
)

speed_cmems_ano = np.sqrt(ugosa.values**2 + vgosa.values**2)


# ============================================================
# 13. ADT + 絶対地衡流 を重ね描き
# ============================================================

overlay_two_vector_fields(
    background=adt.values,
    u1=ugos.values,
    v1=vgos.values,
    u2=ug_adt_main,
    v2=vg_adt_main,
    title=f"ADT and Absolute Geostrophic Velocity\n{date_str}  ({method_main} difference)",
    cbar_label="ADT (m)",
    cmap=ssh_cmap,
    contour_levels=np.arange(0.2, 2.1, 0.1),
    vmin=adt_vmin,
    vmax=adt_vmax
)


# ============================================================
# 14. SLA + 地衡流偏差 を重ね描き
# ============================================================

overlay_two_vector_fields(
    background=sla.values,
    u1=ugosa.values,
    v1=vgosa.values,
    u2=ug_sla_main,
    v2=vg_sla_main,
    title=f"SLA and Geostrophic Velocity Anomaly\n{date_str}  ({method_main} difference)",
    cbar_label="SLA (m)",
    cmap=sla_cmap,
    contour_levels=np.arange(-0.5, 0.55, 0.05),
    norm=TwoSlopeNorm(vmin=-sla_abs, vcenter=0.0, vmax=sla_abs)
)


# ============================================================
# 15. CMEMSとの比較統計
# ============================================================

du_abs, dv_abs, dspeed_abs = compare_velocity(
    ug_adt_main, vg_adt_main,
    ugos.values, vgos.values,
    f"ADT-derived velocity ({method_main}) vs CMEMS ugos/vgos"
)

du_ano, dv_ano, dspeed_ano = compare_velocity(
    ug_sla_main, vg_sla_main,
    ugosa.values, vgosa.values,
    f"SLA-derived velocity anomaly ({method_main}) vs CMEMS ugosa/vgosa"
)


# ============================================================
# 16. 差分マップ
# ============================================================

# ADT由来の絶対地衡流との差
plot_difference_map(
    du_abs,
    title=f"Difference in zonal velocity\nCalculated from ADT - CMEMS ugos\n{date_str}",
    label="du (m/s)",
    vmax=0.10
)

plot_difference_map(
    dv_abs,
    title=f"Difference in meridional velocity\nCalculated from ADT - CMEMS vgos\n{date_str}",
    label="dv (m/s)",
    vmax=0.10
)

plot_difference_map(
    dspeed_abs,
    title=f"Difference in speed\nCalculated from ADT - CMEMS speed\n{date_str}",
    label="speed difference (m/s)",
    vmax=0.10
)

# SLA由来の地衡流偏差との差
plot_difference_map(
    du_ano,
    title=f"Difference in zonal velocity anomaly\nCalculated from SLA - CMEMS ugosa\n{date_str}",
    label="du anomaly (m/s)",
    vmax=0.02
)

plot_difference_map(
    dv_ano,
    title=f"Difference in meridional velocity anomaly\nCalculated from SLA - CMEMS vgosa\n{date_str}",
    label="dv anomaly (m/s)",
    vmax=0.02
)

plot_difference_map(
    dspeed_ano,
    title=f"Difference in anomaly speed\nCalculated from SLA - CMEMS anomaly speed\n{date_str}",
    label="speed difference (m/s)",
    vmax=0.02
)


# ============================================================
# 17. 中央差分と片側差分の比較（任意）
# ============================================================

if compare_two_methods:

    # ADT
    ug_adt_c, vg_adt_c, speed_adt_c = calc_geostrophic_current_switch(
        adt.values, lon, lat, method="central"
    )
    ug_adt_o, vg_adt_o, speed_adt_o = calc_geostrophic_current_switch(
        adt.values, lon, lat, method="onesided"
    )

    # SLA
    ug_sla_c, vg_sla_c, speed_sla_c = calc_geostrophic_current_switch(
        sla.values, lon, lat, method="central"
    )
    ug_sla_o, vg_sla_o, speed_sla_o = calc_geostrophic_current_switch(
        sla.values, lon, lat, method="onesided"
    )

    # 統計比較
    compare_velocity(
        ug_adt_c, vg_adt_c,
        ugos.values, vgos.values,
        "ADT-derived velocity: central difference vs CMEMS ugos/vgos"
    )

    compare_velocity(
        ug_adt_o, vg_adt_o,
        ugos.values, vgos.values,
        "ADT-derived velocity: one-sided difference vs CMEMS ugos/vgos"
    )

    compare_velocity(
        ug_sla_c, vg_sla_c,
        ugosa.values, vgosa.values,
        "SLA-derived velocity anomaly: central difference vs CMEMS ugosa/vgosa"
    )

    compare_velocity(
        ug_sla_o, vg_sla_o,
        ugosa.values, vgosa.values,
        "SLA-derived velocity anomaly: one-sided difference vs CMEMS ugosa/vgosa"
    )

    # 差分法どうしの差
    du_method_adt = ug_adt_o - ug_adt_c
    dv_method_adt = vg_adt_o - vg_adt_c
    dspeed_method_adt = speed_adt_o - speed_adt_c

    du_method_sla = ug_sla_o - ug_sla_c
    dv_method_sla = vg_sla_o - vg_sla_c
    dspeed_method_sla = speed_sla_o - speed_sla_c

    # ADTの差分法比較
    plot_difference_map(
        du_method_adt,
        title=f"Difference in zonal velocity\none-sided - central (ADT)\n{date_str}",
        label="du difference (m/s)",
        vmax=0.15
    )

    plot_difference_map(
        dv_method_adt,
        title=f"Difference in meridional velocity\none-sided - central (ADT)\n{date_str}",
        label="dv difference (m/s)",
        vmax=0.15
    )

    plot_difference_map(
        dspeed_method_adt,
        title=f"Difference in speed\none-sided - central (ADT)\n{date_str}",
        label="speed difference (m/s)",
        vmax=0.15
    )

    # SLAの差分法比較
    plot_difference_map(
        du_method_sla,
        title=f"Difference in zonal velocity anomaly\none-sided - central (SLA)\n{date_str}",
        label="du anomaly difference (m/s)",
        vmax=0.02
    )

    plot_difference_map(
        dv_method_sla,
        title=f"Difference in meridional velocity anomaly\none-sided - central (SLA)\n{date_str}",
        label="dv anomaly difference (m/s)",
        vmax=0.02
    )

    plot_difference_map(
        dspeed_method_sla,
        title=f"Difference in speed anomaly\none-sided - central (SLA)\n{date_str}",
        label="speed anomaly difference (m/s)",
        vmax=0.02
    )


# ============================================================
# 18. 終了処理
# ============================================================

ds.close()
print("Finished.")

7. 結果の読み方

中央差分を使うと、SLAから計算した地衡流偏差は、CMEMSの ugosa, vgosa と非常によく一致しやすい。一方、ADTから計算した絶対地衡流は、黒潮流軸や沿岸付近でCMEMSの ugos, vgos と差が残りやすい。

比較	期待される特徴	考えること
SLA → ugosa/vgosa	中央差分ではほぼ一致しやすい	地衡流偏差はSLAの勾配から計算されることを確認できる
ADT → ugos/vgos	大局的には一致するが、黒潮流軸・沿岸で差が残りやすい	ADTにはMDTが含まれ、プロダクト処理、平滑化、マスク処理の影響を受ける
片側差分 vs 中央差分	片側差分の方が差が大きくなりやすい	数値微分の方法が流速推定に影響する

考えてみよう

黒潮流軸付近で差が大きくなるのはなぜか。
なぜSLA由来の偏差流速はCMEMSとよく一致しやすいのか。
ADT由来の絶対流速が完全一致しない理由を、MDT、平滑化、マスク処理という言葉を使って説明せよ。
観測実習海域は黒潮流軸の北側か南側か。地図から判断せよ。

8. まとめ

海面高度の水平勾配から地衡流を計算できる。
地衡流計算では、コリオリパラメータ、距離、単位変換、差分法が重要である。
中央差分と片側差分では、黒潮のような強勾配域で結果が変わる。
CMEMSの流速を盲目的に使うのではなく、変数の意味と計算方法を確認する必要がある。