11.5と同じ黒潮海域のSLAデータを使い、dx、dy、SLA勾配、コリオリパラメータ、地衡流、相対渦度、EKEを穴埋め形式で1つずつ計算する。
前のページでは、黒潮海域のSLAから地衡流、相対渦度、EKEを段階的に計算した。このページでは同じ処理を、学生が自分で空欄を埋めながら実行する。
ポイントは、SLA → 距離 dx, dy → 勾配 dη/dx, dη/dy → 地衡流 ug, vg → 流速の微分 → 相対渦度 ζ → EKE という変換の流れを理解することである。
このページでは、11と同じ直近の黒潮海域衛星海面高度データを使う。以前の12で使っていた月平均SLAファイルではなく、11.5と同じ日別データに統一する。
| 変数 | 意味 | このページでの使い方 |
|---|---|---|
sla | Sea Level Anomaly | SLA勾配から偏差地衡流を計算する。 |
ugosa, vgosa | CMEMSに含まれるSLA由来の地衡流偏差 | 自分で計算した ug, vg と比較する。 |
longitude, latitude | 経度・緯度 | dx, dy, f の計算に使う。 |
ugosa, vgosa と比較する。東西流速 ug は南北方向のSLA勾配 dη/dy から、南北流速 vg は東西方向のSLA勾配 dη/dx から計算される。
dyは南北方向の中央差分距離であり、今回の格子ではほぼ一定である。カラーバーが非常に小さな差を拡大するため、縞模様に見えることがあるが、計算がおかしいわけではない。
dyはほぼ一定である。一方、dxは cos(latitude) が入るため、高緯度ほど短くなる。
授業時間が足りない場合は、STEP 10の ug, vg までで止めてもよい。相対渦度とEKEはその発展として扱う。
各STEPの空欄 ____ を埋めながら実行する。すべてを一気に貼るのではなく、Jupyter Labでセルごとに実行して、図を確認する。
# ============================================================
# 12. SLAから地衡流・相対渦度・EKEを段階的に計算する（穴埋め版）
# ============================================================
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 TwoSlopeNorm
# 入力ファイル
filename = "____________________________________.nc"
# 最後の日を使う
it = ____
# 地図範囲
lon_min, lon_max = ____, ____
lat_min, lat_max = ____, ____
# 定数
g = ____
Omega = __________
R = __________
# ベクトルの間引き
skip = ____# ============================================================
# STEP 1. データ読み込み
# ============================================================
ds = xr.open_dataset(________)
print(ds)
lon = ds["_________"].values
lat = ds["________"].values
date_str = str(ds["____"].values[it])[:10]
print("Selected date:", date_str)
# 今回はSLAを使う
eta = ds["____"].isel(time=it).values
# 比較用：CMEMSに入っているSLA由来の地衡流偏差
ugosa = ds["_____"].isel(time=it).values
vgosa = ds["_____"].isel(time=it).values
print("eta shape:", eta.shape)
print("lon shape:", lon.shape)
print("lat shape:", lat.shape)# ============================================================
# STEP 2. 地図描画用関数
# ============================================================
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
def plot_map(field, title, cbar_label, cmap="viridis",
vmin=None, vmax=None, center_zero=False):
fig = plt.figure(figsize=(10, 7))
ax = plt.axes(projection=ccrs.PlateCarree())
setup_map(ax)
if center_zero:
if vmax is None:
vmax = np.nanpercentile(np.abs(field), ____)
norm = TwoSlopeNorm(vmin=-vmax, vcenter=0.0, vmax=vmax)
p = ax.pcolormesh(
lon, lat, field,
transform=ccrs.PlateCarree(),
shading="auto",
cmap=cmap,
norm=norm
)
else:
p = ax.pcolormesh(
lon, lat, field,
transform=ccrs.PlateCarree(),
shading="auto",
cmap=cmap,
vmin=vmin,
vmax=vmax
)
cb = plt.colorbar(p, ax=ax, shrink=0.85, pad=0.03)
cb.set_label(cbar_label)
plt.title(title)
plt.show()# ============================================================
# STEP 3. SLAそのものを可視化
# ============================================================
plot_map(
eta,
title=f"SLA\n{date_str}",
cbar_label="SLA (m)",
cmap="________",
center_zero=____
)# ============================================================
# STEP 4. 緯度経度をラジアンに変換
# ============================================================
lon_rad = np.deg2rad(____)
lat_rad = np.deg2rad(____)
lon2d, lat2d = np.meshgrid(____, ____)
lat2d_rad = np.deg2rad(______)
print("lon range:", lon.min(), lon.max())
print("lat range:", lat.min(), lat.max())
print("lon_rad range:", lon_rad.min(), lon_rad.max())
print("lat_rad range:", lat_rad.min(), lat_rad.max())# ============================================================
# STEP 5. dxを計算する
# ============================================================
# 中央差分では、点 j の右隣 j+1 と左隣 j-1 を使う。
# dx = R cos(phi) [lambda(j+1) - lambda(j-1)]
# ============================================================
dlon2 = lon_rad[____] - lon_rad[____]
dx = np.full_like(eta, np.nan, dtype=float)
dx[:, 1:-1] = (
R
* np.cos(lat2d_rad[:, 1:-1])
* dlon2[np.newaxis, :]
)
print("dx shape:", dx.shape)
print("dx min/max [km]:", np.nanmin(dx)/1000, np.nanmax(dx)/1000)
plot_map(
dx / ____,
title="dx for central difference",
cbar_label="dx (km)",
cmap="viridis"
)# ============================================================
# STEP 6. dyを計算する
# ============================================================
# 中央差分では、点 i の上側 i+1 と下側 i-1 を使う。
# dy = R [phi(i+1) - phi(i-1)]
# ============================================================
dlat2 = lat_rad[____] - lat_rad[____]
dy = np.full_like(eta, np.nan, dtype=float)
dy[1:-1, :] = R * dlat2[:, np.newaxis]
print("dy shape:", dy.shape)
print("dy min/max [km]:", np.nanmin(dy)/1000, np.nanmax(dy)/1000)
plot_map(
dy / ____,
title="dy for central difference",
cbar_label="dy (km)",
cmap="viridis"
)# ============================================================
# STEP 7. 東西方向のSLA勾配 dη/dx を計算する
# ============================================================
# dη/dx = [η(j+1) - η(j-1)] / [x(j+1) - x(j-1)]
# ============================================================
deta_dx = np.full_like(eta, np.nan, dtype=float)
deta_dx[:, 1:-1] = (eta[:, ____] - eta[:, ____]) / dx[:, 1:-1]
print("deta_dx min/max:", np.nanmin(deta_dx), np.nanmax(deta_dx))
plot_map(
deta_dx * ____,
title="SLA gradient in x direction: dη/dx",
cbar_label="dη/dx (×10⁻⁶)",
cmap="RdBu_r",
center_zero=True
)# ============================================================
# STEP 8. 南北方向のSLA勾配 dη/dy を計算する
# ============================================================
# dη/dy = [η(i+1) - η(i-1)] / [y(i+1) - y(i-1)]
# ============================================================
deta_dy = np.full_like(eta, np.nan, dtype=float)
deta_dy[1:-1, :] = (eta[____, :] - eta[____, :]) / dy[1:-1, :]
print("deta_dy min/max:", np.nanmin(deta_dy), np.nanmax(deta_dy))
plot_map(
deta_dy * ____,
title="SLA gradient in y direction: dη/dy",
cbar_label="dη/dy (×10⁻⁶)",
cmap="RdBu_r",
center_zero=True
)# ============================================================
# STEP 9. コリオリパラメータ f を計算する
# ============================================================
# f = 2Ωsinφ
# ============================================================
f = ____ * Omega * np.sin(________)
print("f min/max:", np.nanmin(f), np.nanmax(f))
plot_map(
f * ____,
title="Coriolis parameter f",
cbar_label="f (×10⁻⁴ s⁻¹)",
cmap="viridis"
)# ============================================================
# STEP 10. 地衡流を計算する
# ============================================================
# ug = - (g/f) dη/dy
# vg = (g/f) dη/dx
# ============================================================
ug = - (g / f) * ________
vg = (g / f) * ________
speed = np.sqrt(ug**2 + vg**2)
print("ug min/max:", np.nanmin(ug), np.nanmax(ug))
print("vg min/max:", np.nanmin(vg), np.nanmax(vg))
print("speed min/max:", np.nanmin(speed), np.nanmax(speed))
plot_map(
ug,
title="SLA-derived zonal geostrophic velocity ug",
cbar_label="ug (m/s)",
cmap="RdBu_r",
center_zero=True
)
plot_map(
vg,
title="SLA-derived meridional geostrophic velocity vg",
cbar_label="vg (m/s)",
cmap="RdBu_r",
center_zero=True
)
plot_map(
speed,
title="SLA-derived geostrophic speed",
cbar_label="speed (m/s)",
cmap="turbo",
vmin=0,
vmax=np.nanpercentile(speed, 98)
)# ============================================================
# STEP 11. SLAと地衡流ベクトルを重ねる
# ============================================================
fig = plt.figure(figsize=(10, 7))
ax = plt.axes(projection=ccrs.PlateCarree())
setup_map(ax)
sla_abs = np.nanpercentile(np.abs(eta), ____)
p = ax.pcolormesh(
lon, lat, eta,
transform=ccrs.PlateCarree(),
shading="auto",
cmap="RdBu_r",
norm=TwoSlopeNorm(vmin=-sla_abs, vcenter=0.0, vmax=sla_abs)
)
ax.contour(
lon, lat, eta,
levels=np.arange(-0.5, 0.55, 0.05),
colors="black",
linewidths=0.45,
alpha=0.6,
transform=ccrs.PlateCarree()
)
q = ax.quiver(
lon[::skip],
lat[::skip],
ug[::skip, ::skip],
vg[::skip, ::skip],
transform=ccrs.PlateCarree(),
color="black",
scale=12,
width=0.002
)
ax.quiverkey(
q,
X=0.85,
Y=-0.08,
U=____,
label="0.5 m/s",
labelpos="E",
coordinates="axes"
)
cb = plt.colorbar(p, ax=ax, shrink=0.85, pad=0.03)
cb.set_label("SLA (m)")
plt.title(f"SLA and SLA-derived geostrophic velocity\n{date_str}")
plt.show()# ============================================================
# STEP 12. CMEMSの地衡流偏差と比較する
# ============================================================
du = ____ - ugosa
dv = ____ - vgosa
plot_map(
du,
title="Difference in ug: calculated - CMEMS ugosa",
cbar_label="du (m/s)",
cmap="RdBu_r",
center_zero=True,
vmax=0.02
)
plot_map(
dv,
title="Difference in vg: calculated - CMEMS vgosa",
cbar_label="dv (m/s)",
cmap="RdBu_r",
center_zero=True,
vmax=0.02
)# ============================================================
# STEP 13. ∂v/∂x を計算する
# ============================================================
dvdx = np.full_like(vg, np.nan, dtype=float)
dvdx[:, 1:-1] = (vg[:, ____] - vg[:, ____]) / dx[:, 1:-1]
print("dvdx min/max:", np.nanmin(dvdx), np.nanmax(dvdx))
plot_map(
dvdx * ____,
title="∂v/∂x",
cbar_label="∂v/∂x (×10⁻⁶ s⁻¹)",
cmap="RdBu_r",
center_zero=True
)# ============================================================
# STEP 14. ∂u/∂y を計算する
# ============================================================
dudy = np.full_like(ug, np.nan, dtype=float)
dudy[1:-1, :] = (ug[____, :] - ug[____, :]) / dy[1:-1, :]
print("dudy min/max:", np.nanmin(dudy), np.nanmax(dudy))
plot_map(
dudy * ____,
title="∂u/∂y",
cbar_label="∂u/∂y (×10⁻⁶ s⁻¹)",
cmap="RdBu_r",
center_zero=True
)# ============================================================
# STEP 15. 相対渦度 ζ と ζ/f を計算する
# ============================================================
# ζ = ∂v/∂x - ∂u/∂y
# ============================================================
zeta = ______ - ______
plot_map(
zeta * ____,
title="Relative vorticity ζ",
cbar_label="ζ (×10⁻⁶ s⁻¹)",
cmap="RdBu_r",
center_zero=True
)
zeta_over_f = zeta / ____
plot_map(
zeta_over_f,
title="Relative vorticity normalized by f: ζ/f",
cbar_label="ζ/f",
cmap="RdBu_r",
center_zero=True,
vmax=0.5
)# ============================================================
# STEP 16. EKEを計算する
# ============================================================
# EKE = 1/2 (u^2 + v^2)
# ============================================================
eke = ____ * (ug**2 + vg**2)
plot_map(
eke,
title="Eddy kinetic energy from SLA-derived geostrophic velocity",
cbar_label="EKE (m²/s²)",
cmap="turbo",
vmin=0,
vmax=np.nanpercentile(eke, 98)
)
log10_eke = np.log10(eke)
plot_map(
log10_eke,
title="log10(EKE)",
cbar_label="log10(EKE) [m²/s²]",
cmap="turbo",
vmin=np.nanpercentile(log10_eke, 2),
vmax=np.nanpercentile(log10_eke, 98)
)
ds.close()
print("Finished.")ugosa, vgosa と完全一致しない理由を考えよ。課題提出後、確認用としてパスワードを入力すると、穴埋めの要点と完全スクリプトを表示する。
kuroshio_altimetry_daily_20260501_20260512
-1
100, 150
20, 50
9.81
7.2921e-5
6371000.0
16
filename
longitude
latitude
time
sla
ugosa
vgosa
98
RdBu_r
True
lon, lat
lon, lat
lat2d
2:, :-2
1000
2:, :-2
1000
2:, :-2
1e6
2:, :-2
1e6
2.0, lat2d_rad
1e4
deta_dy
deta_dx
98
0.5
ug
vg
2:, :-2
1e6
2:, :-2
1e6
dvdx, dudy
1e6
f
0.5# ============================================================
# 12. SLAから地衡流・相対渦度・EKEを段階的に計算する（完全版）
# ============================================================
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 TwoSlopeNorm
filename = "kuroshio_altimetry_daily_20260501_20260512.nc"
it = -1
lon_min, lon_max = 100, 150
lat_min, lat_max = 20, 50
g = 9.81
Omega = 7.2921e-5
R = 6371000.0
skip = 16
# -----------------------------
# Read data
# -----------------------------
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)
eta = ds["sla"].isel(time=it).values
ugosa = ds["ugosa"].isel(time=it).values
vgosa = ds["vgosa"].isel(time=it).values
print("eta shape:", eta.shape)
print("lon shape:", lon.shape)
print("lat shape:", lat.shape)
# -----------------------------
# Map helpers
# -----------------------------
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
def plot_map(field, title, cbar_label, cmap="viridis", vmin=None, vmax=None, center_zero=False):
fig = plt.figure(figsize=(10, 7))
ax = plt.axes(projection=ccrs.PlateCarree())
setup_map(ax)
if center_zero:
if vmax is None:
vmax = np.nanpercentile(np.abs(field), 98)
norm = TwoSlopeNorm(vmin=-vmax, vcenter=0.0, vmax=vmax)
p = ax.pcolormesh(lon, lat, field, transform=ccrs.PlateCarree(), shading="auto", cmap=cmap, norm=norm)
else:
p = ax.pcolormesh(lon, lat, field, transform=ccrs.PlateCarree(), shading="auto", cmap=cmap, vmin=vmin, vmax=vmax)
cb = plt.colorbar(p, ax=ax, shrink=0.85, pad=0.03)
cb.set_label(cbar_label)
plt.title(title)
plt.show()
# -----------------------------
# SLA
# -----------------------------
plot_map(eta, title=f"SLA\n{date_str}", cbar_label="SLA (m)", cmap="RdBu_r", center_zero=True)
# -----------------------------
# Coordinates and metric
# -----------------------------
lon_rad = np.deg2rad(lon)
lat_rad = np.deg2rad(lat)
lon2d, lat2d = np.meshgrid(lon, lat)
lat2d_rad = np.deg2rad(lat2d)
# dx
dlon2 = lon_rad[2:] - lon_rad[:-2]
dx = np.full_like(eta, np.nan, dtype=float)
dx[:, 1:-1] = R * np.cos(lat2d_rad[:, 1:-1]) * dlon2[np.newaxis, :]
print("dx shape:", dx.shape)
print("dx min/max [km]:", np.nanmin(dx)/1000, np.nanmax(dx)/1000)
plot_map(dx / 1000, title="dx for central difference", cbar_label="dx (km)", cmap="viridis")
# dy
dlat2 = lat_rad[2:] - lat_rad[:-2]
dy = np.full_like(eta, np.nan, dtype=float)
dy[1:-1, :] = R * dlat2[:, np.newaxis]
print("dy shape:", dy.shape)
print("dy min/max [km]:", np.nanmin(dy)/1000, np.nanmax(dy)/1000)
plot_map(dy / 1000, title="dy for central difference", cbar_label="dy (km)", cmap="viridis")
# -----------------------------
# SLA gradients
# -----------------------------
deta_dx = np.full_like(eta, np.nan, dtype=float)
deta_dx[:, 1:-1] = (eta[:, 2:] - eta[:, :-2]) / dx[:, 1:-1]
plot_map(deta_dx * 1e6, title="SLA gradient in x direction: dη/dx", cbar_label="dη/dx (×10⁻⁶)", cmap="RdBu_r", center_zero=True)
deta_dy = np.full_like(eta, np.nan, dtype=float)
deta_dy[1:-1, :] = (eta[2:, :] - eta[:-2, :]) / dy[1:-1, :]
plot_map(deta_dy * 1e6, title="SLA gradient in y direction: dη/dy", cbar_label="dη/dy (×10⁻⁶)", cmap="RdBu_r", center_zero=True)
# -----------------------------
# Coriolis parameter
# -----------------------------
f = 2.0 * Omega * np.sin(lat2d_rad)
plot_map(f * 1e4, title="Coriolis parameter f", cbar_label="f (×10⁻⁴ s⁻¹)", cmap="viridis")
# -----------------------------
# Geostrophic velocity
# -----------------------------
ug = - (g / f) * deta_dy
vg = (g / f) * deta_dx
speed = np.sqrt(ug**2 + vg**2)
plot_map(ug, title="SLA-derived zonal geostrophic velocity ug", cbar_label="ug (m/s)", cmap="RdBu_r", center_zero=True)
plot_map(vg, title="SLA-derived meridional geostrophic velocity vg", cbar_label="vg (m/s)", cmap="RdBu_r", center_zero=True)
plot_map(speed, title="SLA-derived geostrophic speed", cbar_label="speed (m/s)", cmap="turbo", vmin=0, vmax=np.nanpercentile(speed, 98))
# SLA + vectors
fig = plt.figure(figsize=(10, 7))
ax = plt.axes(projection=ccrs.PlateCarree())
setup_map(ax)
sla_abs = np.nanpercentile(np.abs(eta), 98)
p = ax.pcolormesh(lon, lat, eta, transform=ccrs.PlateCarree(), shading="auto", cmap="RdBu_r", norm=TwoSlopeNorm(vmin=-sla_abs, vcenter=0.0, vmax=sla_abs))
ax.contour(lon, lat, eta, levels=np.arange(-0.5, 0.55, 0.05), colors="black", linewidths=0.45, alpha=0.6, transform=ccrs.PlateCarree())
q = ax.quiver(lon[::skip], lat[::skip], ug[::skip, ::skip], vg[::skip, ::skip], transform=ccrs.PlateCarree(), color="black", scale=12, width=0.002)
ax.quiverkey(q, X=0.85, Y=-0.08, U=0.5, label="0.5 m/s", labelpos="E", coordinates="axes")
cb = plt.colorbar(p, ax=ax, shrink=0.85, pad=0.03)
cb.set_label("SLA (m)")
plt.title(f"SLA and SLA-derived geostrophic velocity\n{date_str}")
plt.show()
# Difference from CMEMS anomaly velocity
du = ug - ugosa
dv = vg - vgosa
plot_map(du, title="Difference in ug: calculated - CMEMS ugosa", cbar_label="du (m/s)", cmap="RdBu_r", center_zero=True, vmax=0.02)
plot_map(dv, title="Difference in vg: calculated - CMEMS vgosa", cbar_label="dv (m/s)", cmap="RdBu_r", center_zero=True, vmax=0.02)
# -----------------------------
# Relative vorticity
# -----------------------------
dvdx = np.full_like(vg, np.nan, dtype=float)
dvdx[:, 1:-1] = (vg[:, 2:] - vg[:, :-2]) / dx[:, 1:-1]
plot_map(dvdx * 1e6, title="∂v/∂x", cbar_label="∂v/∂x (×10⁻⁶ s⁻¹)", cmap="RdBu_r", center_zero=True)
dudy = np.full_like(ug, np.nan, dtype=float)
dudy[1:-1, :] = (ug[2:, :] - ug[:-2, :]) / dy[1:-1, :]
plot_map(dudy * 1e6, title="∂u/∂y", cbar_label="∂u/∂y (×10⁻⁶ s⁻¹)", cmap="RdBu_r", center_zero=True)
zeta = dvdx - dudy
plot_map(zeta * 1e6, title="Relative vorticity ζ", cbar_label="ζ (×10⁻⁶ s⁻¹)", cmap="RdBu_r", center_zero=True)
zeta_over_f = zeta / f
plot_map(zeta_over_f, title="Relative vorticity normalized by f: ζ/f", cbar_label="ζ/f", cmap="RdBu_r", center_zero=True, vmax=0.5)
# -----------------------------
# EKE
# -----------------------------
eke = 0.5 * (ug**2 + vg**2)
plot_map(eke, title="Eddy kinetic energy from SLA-derived geostrophic velocity", cbar_label="EKE (m²/s²)", cmap="turbo", vmin=0, vmax=np.nanpercentile(eke, 98))
log10_eke = np.log10(eke)
plot_map(log10_eke, title="log10(EKE)", cbar_label="log10(EKE) [m²/s²]", cmap="turbo", vmin=np.nanpercentile(log10_eke, 2), vmax=np.nanpercentile(log10_eke, 98))
ds.close()
print("Finished.")