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

4.1 まず「微分」を格子点の「差」で考える

地衡流の式には ∂η/∂x や ∂η/∂y が出てくる。これは、海面高度 η が東西方向・南北方向にどれくらい傾いているかを表す。

∂η/∂x ≃ 東西方向の海面高度差 / 東西方向の距離

∂η/∂y ≃ 南北方向の海面高度差 / 南北方向の距離

しかし、衛星データは連続した曲線ではなく、緯度経度の格子点上に並んだ値である。そのため、微分を直接計算するのではなく、隣り合う格子点の差で近似する。

η0   η1   η2   η3   η4
●----●----●----●----●
x0   x1   x2   x3   x4

たとえば、点 x2 付近の傾きを知りたい場合、中央差分では、点 x2 をはさんだ左右の値を使う。

点 x2 の傾き ≃ (η3 − η1) / (x3 − x1)

中央差分では、注目点そのものではなく、注目点をはさんだ左右または上下の値を使って傾きを近似する。

4.2 Pythonのスライス 2:, :-2, 1:-1 の意味

中央差分のコードでわかりにくいのは、次のような配列の切り出しである。

dlon2 = lon_rad[2:] - lon_rad[:-2]
d_eta_dx[:, 1:-1] = (eta[:, 2:] - eta[:, :-2]) / dx2d

まず、1次元の配列で考える。

eta = [10, 12, 15, 21, 22, 25]

書き方	取り出される値	意味
eta[2:]	[15, 21, 22, 25]	2番目以降。中央差分では「右隣側」の値に対応する。
eta[:-2]	[10, 12, 15, 21]	最後の2個を除く。中央差分では「左隣側」の値に対応する。
eta[1:-1]	[12, 15, 21, 22]	最初と最後を除いた内側。計算結果を入れる中央の点に対応する。

したがって、

eta[2:] - eta[:-2]

は、次の差をまとめて計算している。

[15-10, 21-12, 22-15, 25-21]

これは、各点について η[i+1] − η[i−1] を計算していることになる。

中央の点	使う差	意味
i = 1	η[2] − η[0]	点1をはさんだ右と左の差
i = 2	η[3] − η[1]	点2をはさんだ右と左の差
i = 3	η[4] − η[2]	点3をはさんだ右と左の差
i = 4	η[5] − η[3]	点4をはさんだ右と左の差

覚え方：2: は右隣側、:-2 は左隣側、1:-1 は結果を入れる中央の場所、と考えるとよい。

4.3 なぜ最初と最後の格子点は計算できないのか

中央差分では、注目点の左隣と右隣が必要である。したがって、最初の点には左隣がなく、最後の点には右隣がない。

0   1   2   3   4   5
●---●---●---●---●---●
    ↑   ↑   ↑   ↑
  計算できる点

そのため、中央差分の結果は、元の配列と同じ大きさの配列を NaN で用意しておき、内側だけに値を入れている。

d_eta_dx = np.full_like(eta, np.nan, dtype=float)
d_eta_dx[:, 1:-1] = (eta[:, 2:] - eta[:, :-2]) / dx2d

1:-1 は「最初と最後を除いた内側」を意味する。外側1格子は中央差分では計算できないので、NaN のまま残る。

4.4 2次元配列では、東西方向と南北方向でスライスが変わる

海面高度 eta は、緯度 × 経度の2次元配列である。一般に、

eta[緯度方向, 経度方向]

のように考える。

計算したい勾配	Pythonの書き方	何をしているか
東西方向の勾配 dη/dx	eta[:, 2:] - eta[:, :-2]	すべての緯度について、経度方向に右隣 − 左隣を計算する。
南北方向の勾配 dη/dy	eta[2:, :] - eta[:-2, :]	すべての経度について、緯度方向に上側 − 下側を計算する。

つまり、コロン : が入っている方向は「全部使う」という意味である。

# x方向、つまり経度方向の差
eta[:, 2:] - eta[:, :-2]

# y方向、つまり緯度方向の差
eta[2:, :] - eta[:-2, :]

2次元配列では、どちらの軸をずらしているかを意識することが重要である。

4.5 緯度経度の差をメートルに直す

NetCDFの座標は、たとえば longitude = 139.0, 139.25, 139.5 のように「度」で表されている。しかし地衡流の式に出てくる x, y は、メートル単位の距離である。

緯度経度のまま差を取ると、海面高度勾配の単位は m/degree になってしまう。地衡流計算に必要なのは m/m の勾配である。

そのため、緯度経度の差をラジアンに直し、地球半径をかけて距離に変換する。

lon_rad = np.deg2rad(lon)
lat_rad = np.deg2rad(lat)

南北方向の距離は、地球半径を R とすると、

dy = R dφ

で計算できる。ここで dφ は緯度差をラジアンで表したものである。

dlat2 = lat_rad[2:] - lat_rad[:-2]
dy2d = R * dlat2[:, np.newaxis]

一方、東西方向の距離は緯度によって変わる。赤道付近では経度1度の距離は長く、高緯度では短くなる。そのため、緯度の余弦 cos(latitude) をかける。

dx = R cos(φ) dλ

dlon2 = lon_rad[2:] - lon_rad[:-2]
dx2d = R * np.cos(lat2d_rad[:, 1:-1]) * dlon2[np.newaxis, :]

方向	距離計算	理由
南北方向	dy = R dlat	緯度1度あたりの距離はほぼ一定。
東西方向	dx = R cos(lat) dlon	経度1度あたりの距離は緯度によって変わる。

4.6 np.newaxis は何をしているか

np.newaxis は、配列の次元を1つ増やすための書き方である。ここでは、1次元の距離配列を、2次元の緯度経度格子と掛け算できる形にするために使っている。

dlon2.shape
# 例: (198,)

dlon2[np.newaxis, :].shape
# 例: (1, 198)

dlon2 は経度方向だけの1次元配列である。一方、np.cos(lat2d_rad[:, 1:-1]) は緯度 × 経度の2次元配列である。両者を掛け算するために、dlon2[np.newaxis, :] として横方向の2次元配列にしている。

ここでやっていることは、難しい数学ではなく、各緯度・各経度格子に対応する東西距離 dx を作るという処理である。

4.7 このページの中央差分コードを日本語に直す

最終的に、中央差分のコードは次のように読める。

コード	日本語での意味
dlon2 = lon_rad[2:] - lon_rad[:-2]	中央差分で使う、右隣と左隣の経度差を作る。
dx2d = R * cos(lat) * dlon2	経度差を、各緯度での東西距離 m に直す。
eta[:, 2:] - eta[:, :-2]	海面高度について、右隣 − 左隣を計算する。
d_eta_dx[:, 1:-1] = ...	計算した東西勾配を、内側の格子点に入れる。
dlat2 = lat_rad[2:] - lat_rad[:-2]	中央差分で使う、上側と下側の緯度差を作る。
dy2d = R * dlat2[:, np.newaxis]	緯度差を南北距離 m に直し、2次元配列と計算できる形にする。
eta[2:, :] - eta[:-2, :]	海面高度について、上側 − 下側を計算する。
d_eta_dy[1:-1, :] = ...	計算した南北勾配を、内側の格子点に入れる。

この部分の核心：中央差分では、注目点の右隣 − 左隣を、右隣までの距離 − 左隣までの距離で割り、その結果を注目点に入れている。

4.8 片側差分との違い

片側差分では、注目点とその隣の点だけを使う。

片側差分：点 j の傾き ≃ (η[j+1] − η[j]) / (x[j+1] − x[j])

コードでは、次のようになる。

dlon1 = lon_rad[1:] - lon_rad[:-1]
d_eta_dx[:, :-1] = (eta[:, 1:] - eta[:, :-1]) / dx2d

中央差分と比較すると、使っている点が違う。

方法	使う値	特徴
片側差分	η[j+1] − η[j]	端の近くまで計算しやすいが、勾配の位置が格子点の中間にずれやすい。
中央差分	η[j+1] − η[j−1]	注目点を中心に左右対称に見るので、一般に滑らかで精度がよい。ただし端は計算できない。

この授業では、片側差分を「中央差分との違いを見るための比較用」として使う。実際のプロダクトでどの差分・平滑化・マスク処理が使われているかは、データ提供元の仕様を確認する必要がある。

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 の右隣 j+1 と左隣 j-1 を使う。
        # eta[:, 2:]   : 経度方向の右隣側の値
        # eta[:, :-2]  : 経度方向の左隣側の値
        # [:, 1:-1]    : 計算結果を入れる中央の格子点
        # [η(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
        # eta[2:, :]   : 緯度方向の上側の値
        # eta[:-2, :]  : 緯度方向の下側の値
        # [1:-1, :]    : 計算結果を入れる中央の格子点
        # [η(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の流速を盲目的に使うのではなく、変数の意味と計算方法を確認する必要がある。