1. 从零开始：为什么我们需要绘制海温异常图？

如果你关注过天气预报，或者刷到过关于“厄尔尼诺”、“拉尼娜”的新闻，那你可能已经听说过“海温异常”这个词。简单来说，海温异常就是某个时期的海水表面温度，和过去几十年平均水平的差值。听起来好像只是个简单的减法，对吧？但就是这个差值，藏着气候变化的“密码”。

我刚开始接触气象数据可视化的时候，也觉得这玩意儿挺玄乎。一堆枯燥的数字，怎么变成一张能讲故事、能揭示规律的图？后来在项目中实际用起来才发现，一张清晰、专业的海温异常分布图，比任何长篇大论的文字报告都更有说服力。它能直观地告诉你：今年冬天，哪片海洋热得反常，哪片又冷得异常？这些异常区域和我们的台风路径、降雨带、甚至冬天的寒潮有没有关系？

举个例子，2023年冬季，太平洋中东部那片醒目的红色（代表偏暖）区域，就是一次典型的厄尔尼诺事件的“签名”。通过Python，我们可以亲手从庞大的全球数据集中，把这种信号提取并可视化出来。这个过程不仅是对数据分析和编程能力的锻炼，更是理解我们脚下这个星球气候脉搏的绝佳方式。所以，无论你是气象、海洋专业的学生，还是对数据科学和地球科学感兴趣的开发者，掌握这项技能都绝对“入股不亏”。接下来，我就手把手带你，用Python中最强大的地理绘图“黄金搭档”——Cartopy和Matplotlib，从下载数据开始，一步步绘制出一张属于自己的、专业的全球海温异常分布图。

2. 环境搭建与数据准备：磨刀不误砍柴工

工欲善其事，必先利其器。在开始写代码画图之前，我们需要先把“厨房”收拾好，把“食材”准备好。这里说的“厨房”就是你的Python环境，“食材”就是我们要用的海温数据。

2.1 安装必备的Python库

首先，确保你的Python环境（我强烈推荐使用Anaconda来管理，能省去无数依赖冲突的麻烦）里已经安装了以下几个核心库。打开你的终端或命令提示符，一行命令就能搞定：

pip install numpy matplotlib xarray cartopy netCDF4

我来简单介绍一下这几个“得力干将”：

• NumPy：Python科学计算的基石，处理多维数组（我们的海温数据就是一个三维数组：时间 x 纬度 x 经度）离不开它。
• Matplotlib：Python绘图的“老大哥”，功能极其强大，我们用它来创建画布、绘制填色图和颜色条。
• Xarray：这是处理像气象、海洋这种带标签的多维网格数据的“神器”！它比直接用NumPy友好太多了，可以像操作字典一样，通过时间、经纬度来切片数据，避免了繁琐的索引计算。我实测下来，用xarray处理NetCDF格式的气象数据，代码能简洁一半以上。
• Cartopy：本次的“主角”之一。它是一个用于地理空间数据可视化的库，可以轻松处理地图投影、添加海岸线、国界等地理要素。没有它，我们画出来的就只是一张没有地理意义的“温度贴图”。
• netCDF4：一个底层库，用于读取我们数据文件的格式（NetCDF）。通常xarray会依赖它。

安装Cartopy时如果遇到问题（特别是在Windows上），可能需要先安装一些地理信息库的依赖，比如GEOS、Proj。最省心的办法就是去Anaconda的官网搜索cartopy，用conda install -c conda-forge cartopy命令安装，conda会帮你自动处理好所有依赖。

2.2 获取ERA5再分析数据

我们的“食材”来自欧洲中期天气预报中心（ECMWF）的ERA5再分析数据集。你可以把它理解为一个用最先进的数值天气预报模型和全球观测数据“回放”出来的地球气候“纪录片”，它提供了自1940年以来，包括海温、气压、风场等在内的大量高精度、高时空一致性的数据。

数据获取步骤：

1. 访问平台：你需要到ECMWF的Climate Data Store (CDS) 官网注册一个免费账户。
2. 搜索数据：在搜索框中输入“ERA5 monthly averaged data on single levels”，选择对应的数据集。
3. 定制数据：在下载界面，你需要选择：
   • 变量：选择“Sea surface temperature”。
   • 时间范围：选择从1940年1月到2023年12月（或你能获取到的最新月份）。
   • 地理区域：选择全球。
   • 格式：选择NetCDF格式。
4. 提交请求：点击下载后，由于数据量较大，CDS会先在后台准备数据，然后通过邮件通知你下载链接。这个过程可能需要等待几分钟到几小时。

我个人的经验是，可以先下载一个时间范围较小的数据（比如仅包含1991-2023年）用于测试代码，等代码调试无误后，再下载完整的数据集进行最终分析。这样能节省大量等待时间。下载完成后，你会得到一个类似era5_monthly_sst_1940_2023.nc的文件，这就是我们后续操作的起点。

3. 核心代码实战：一步步解析数据处理与绘图

数据到手，环境就绪，现在让我们进入最核心的代码环节。我会把每一段代码掰开揉碎了讲，确保你不仅能复制粘贴跑通，更能理解背后的逻辑。

3.1 导入库与读取数据

打开你的Python编辑器（Jupyter Notebook, VS Code, PyCharm都可以），新建一个文件，我们开始写代码。

# 1. 导入必不可少的工具库
import numpy as np
import matplotlib.pyplot as plt
import xarray as xr
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.mpl.ticker as cticker

第一行到第三行是标准的数据处理三件套。第四行ccrs是Cartopy中定义地图投影的模块，比如我们等下要用的PlateCarree（等经纬度投影）。第五行cfeature用于添加地理特征，比如海岸线、陆地、河流。第六行cticker则专门用于格式化地图上的经纬度刻度标签，让它显示为“120°E”而不是单纯的数字120。

# 2. 读取NetCDF数据文件
file_path = './data/ERA5_sst_1940_202307.nc' # 请替换为你的实际文件路径
ds = xr.open_dataset(file_path)
print(ds)

用xr.open_dataset打开文件后，我习惯立刻用print看一下数据集的整体情况。这会输出数据的维度、坐标（经度、纬度、时间）、变量名等信息。你会看到类似这样的结构：

Dimensions:    (time: 1003, latitude: 181, longitude: 360)
Coordinates:
  * time       (time) datetime64[ns] 1940-01-01 1940-02-01 ... 2023-07-01
  * latitude   (latitude) float64 90.0 89.0 88.0 ... -89.0 -90.0
  * longitude  (longitude) float64 0.0 1.0 2.0 ... 357.0 358.0 359.0
Data variables:
    sst        (time, latitude, longitude) float32 ...

看，xarray的优势立刻体现出来了！它清晰地告诉我们，数据变量sst（海表温度）是一个三维数组，维度顺序是(时间， 纬度， 经度)。经纬度是坐标，我们可以直接用名字来索引。

3.2 计算冬季平均与气候态异常

这是数据处理的关键步骤，目的是得到“2023年冬季平均海温”减去“30年冬季气候态平均”的结果。

# 3. 提取冬季月份数据
# 假设我们的时间坐标名为'time'，海温变量名为'sst'
# 提取1991-2020年所有12月、1月、2月的数据
winter_months = [12, 1, 2]
sst_clim = ds.sst.sel(time=ds.time.dt.month.isin(winter_months))
sst_clim = sst_clim.sel(time=slice('1991-01-01', '2020-12-31'))

# 提取2023年冬季数据（即2022年12月，2023年1月、2月）
sst_2023_winter = ds.sst.sel(time=ds.time.dt.month.isin(winter_months))
sst_2023_winter = sst_2023_winter.sel(time=slice('2022-12-01', '2023-02-28'))

# 打印一下看看提取的数据形状，确认无误
print(f"气候态冬季数据形状: {sst_clim.shape}")
print(f"2023年冬季数据形状: {sst_2023_winter.shape}")

这里用了xarray非常强大的.sel()和.isel()方法进行数据选择。ds.time.dt.month.isin(winter_months)这个操作，可以直接筛选出时间维度中月份为12、1、2的所有时间点，非常直观，避免了写循环。slice函数则用于选择时间范围。

# 4. 计算气候态平均和2023年冬季平均
# 沿着时间维度计算平均值
sst_clim_mean = sst_clim.mean(dim='time')
sst_2023_mean = sst_2023_winter.mean(dim='time')

# 5. 计算海温异常
sst_anomaly = sst_2023_mean - sst_clim_mean

# 检查一下结果
print(f"气候态平均SST范围: {sst_clim_mean.min().values:.2f} 到 {sst_clim_mean.max().values:.2f}")
print(f"海温异常范围: {sst_anomaly.min().values:.2f} 到 {sst_anomaly.max().values:.2f}")

计算平均值时，我们指定dim='time'，意思是沿着时间轴求平均。对于sst_clim（假设有30年*3个月=90个月），求平均后得到一个二维数组（纬度 x 经度），这就是30年冬季平均的气候态。对于sst_2023_winter（3个月），求平均后得到2023年冬季的平均海温。两者相减，就得到了我们想要的海温异常。这个异常值如果是正的，表示比常年偏暖；负的则表示偏冷。

3.3 构建地图与绘制填色图

数据处理完毕，现在进入“画图”阶段，这是让数据“活”起来的一步。

# 6. 创建图形和地图投影
fig = plt.figure(figsize=(16, 8)) # 设置一个宽屏的图形大小，适合全球地图

# 关键一步：创建地图投影。central_longitude=180意味着将地图的中心经线设在180度（国际日期变更线附近）
# 这样太平洋就在图中央，更符合我们看全球海温分布的习惯。
proj = ccrs.PlateCarree(central_longitude=180)

# 在图形上添加一个子图，并指定其使用我们定义的地图投影
ax = fig.add_subplot(1, 1, 1, projection=proj)

# 设置地图的显示范围（全球）
ax.set_global()
# 或者更精确地控制：ax.set_extent([-180, 180, -90, 90], crs=proj)

创建投影PlateCarree(central_longitude=180)是让太平洋居中的关键。add_subplot时传入projection=proj参数，告诉Matplotlib这个坐标轴是地理坐标轴。

# 7. 添加地理要素，让地图更丰富
# 添加高分辨率的海岸线
ax.add_feature(cfeature.COASTLINE.with_scale('50m'), linewidth=0.5)
# 添加陆地填充，颜色设为浅灰色，避免和海洋填色混淆
ax.add_feature(cfeature.LAND, color='lightgray')
# 可以添加湖泊、河流等，但为了图面简洁，这里先不加
# ax.add_feature(cfeature.LAKES, alpha=0.5)
# ax.add_feature(cfeature.RIVERS)

# 8. 设置并格式化经纬度网格线
# 设置网格线间隔
gl = ax.gridlines(crs=proj, draw_labels=True, linewidth=0.5, color='gray', alpha=0.5, linestyle='--')
gl.top_labels = False  # 不显示顶部标签
gl.right_labels = False # 不显示右侧标签
gl.xlabel_style = {'size': 10}
gl.ylabel_style = {'size': 10}

cfeature提供了不同精度的自然地理数据。‘50m’精度足够出版级使用。添加陆地并着色，能立刻让地图有层次感。网格线gridlines的draw_labels=True参数非常方便，能自动在合适的侧面显示经纬度标签。

# 9. 绘制海温异常填色图
# 这是最核心的绘图命令
# levels: 定义色阶的分段。np.arange(-3, 3.1, 0.2)表示从-3度到3度，每隔0.2度画一条等值线并填充。
# cmap: 颜色映射。'RdBu_r'是红蓝渐变色，'_r'表示反转，通常让红色代表暖（正异常），蓝色代表冷（负异常）。
# transform: 必须指定！告诉Cartopy数据所在的坐标系统。我们的数据是经纬度，所以用ccrs.PlateCarree()。
# extend: ‘both’表示颜色条两端有箭头，因为数据可能超出我们设定的levels范围。
levels = np.arange(-3, 3.1, 0.2)
cf = ax.contourf(sst_anomaly.longitude, sst_anomaly.latitude, sst_anomaly,
                 levels=levels, cmap='RdBu_r', transform=ccrs.PlateCarree(), extend='both')

# 10. 添加颜色条
# 将颜色条与刚才的填色图cf关联起来
cbar = plt.colorbar(cf, ax=ax, orientation='horizontal', pad=0.05, shrink=0.8)
cbar.set_label('Sea Surface Temperature Anomaly (°C)', fontsize=12) # 设置颜色条标签
cbar.ax.tick_params(labelsize=10) # 设置颜色条刻度字体大小

ax.contourf是绘制填充等值线图的函数。有几个参数我踩过坑：

1. transform=ccrs.PlateCarree()：这个绝对不能少！它声明了输入数据的坐标系。即使你的地图投影是别的（比如Robinson），只要原始数据是经纬度，这里就填PlateCarree。少了它，图画出来位置全是错的。
2. levels：控制填色的精细度和范围。我通常先粗略画一次，看看数据的最大最小值，再确定一个合理的范围。对于海温异常，±3°C是一个比较常见的范围。
3. cmap='RdBu_r'：这是一个在气象学中非常经典的“冷-暖”色系，中间是白色，完美契合“异常”的概念，零值附近是白色，正负异常一目了然。

3.4 添加标题与保存结果

最后，给我们的成果图加上说明，并保存下来。

# 11. 添加标题
ax.set_title('Global Sea Surface Temperature Anomaly (DJF 2022-2023 vs 1991-2020 Climatology)', fontsize=14, pad=20)

# 12. 调整布局并显示图形
plt.tight_layout()
plt.show()

# 13. 保存高清图片
fig.savefig('global_sst_anomaly_2023_winter.png', dpi=300, bbox_inches='tight')
print("图形已保存为 'global_sst_anomaly_2023_winter.png'")

plt.tight_layout()可以自动调整子图参数，使图形元素不重叠，让图更美观。保存时dpi=300确保印刷或报告使用的清晰度，bbox_inches='tight'能裁剪掉图形周围多余的白边。

4. 进阶技巧与美化：让你的图更专业

如果上面的代码跑通，你已经得到了一张合格的海温异常图。但要让它在报告或论文中脱颖而出，还需要一些“美颜”和“精修”。下面分享几个我实践中总结的进阶技巧。

4.1 自定义色阶与零值对齐

Matplotlib内置的‘RdBu_r’色系很好，但有时我们希望更精确地控制，特别是确保0值绝对对应白色。我们可以从‘RdBu_r’色图中截取一部分，或者创建自定义的离散色板。

# 方法一：从连续色图中截取并离散化
import matplotlib.colors as mcolors

# 定义我们想要的色阶范围
level_boundaries = np.arange(-3, 3.1, 0.5)
n_colors = len(level_boundaries) - 1

# 从'RdBu_r'色图中取出对应数量的颜色
# 关键：确保色图的中间颜色（对应0值）被取到。
# cmap(np.linspace(0, 1, n_colors)) 会均匀取色，但0值不一定在正中间。
# 我们需要计算0值在level_boundaries中的相对位置。
zero_index = np.where(level_boundaries <= 0)[0][-1] # 找到0值所在区间的左边界索引
# 构建一个归一化的列表，确保0值两侧的区间数量大致相等（如果区间对称的话）
norm_vals = np.linspace(0, 1, n_colors+1)
custom_cmap_list = plt.cm.RdBu_r(norm_vals)
# 根据0值位置调整，使0值区间对应色图最中间的颜色（索引为len(custom_cmap_list)//2）
# 这里逻辑稍复杂，一个更稳妥的方法是直接使用DivergingNorm（已弃用）或TwoSlopeNorm

更简单有效的方法是使用TwoSlopeNorm进行标准化，它能完美解决“零值居中”问题：

from matplotlib.colors import TwoSlopeNorm

# 定义数据的中心点（vcenter）为0
norm = TwoSlopeNorm(vmin=-3, vcenter=0, vmax=3)
# 绘图时使用这个norm
cf = ax.contourf(sst_anomaly.longitude, sst_anomaly.latitude, sst_anomaly,
                 levels=levels, cmap='RdBu_r', norm=norm, transform=ccrs.PlateCarree(), extend='both')

用了TwoSlopeNorm之后，无论你的levels如何设置，颜色映射都会严格保证vcenter（这里为0）对应色图的中间颜色（‘RdBu_r’的白色部分），正负异常的颜色过渡就非常对称和准确。

4.2 添加显著性检验打点

在科研绘图中，我们不仅关心异常的大小，还关心它是否“显著”，即是否超出了自然变率的范围。通常我们会进行统计检验（如t检验），将通过显著性检验的区域用打点（stippling）或画斜线（hatching）的方式标记出来。

# 假设我们已经计算了气候态30年冬季数据的标准差sst_clim_std（也是一个二维场）
# 并进行了逐格点的t检验，得到了一个布尔型数组`significant`（True表示该格点异常显著）

# 计算标准差
sst_clim_std = sst_clim.std(dim='time')
# 简化示例：这里用一个简单的阈值法模拟显著性区域（实际应用需用t检验）
# 例如，认为异常绝对值超过1.5倍标准差的区域为“显著”
significant = np.abs(sst_anomaly) > (1.5 * sst_clim_std)

# 将显著的格点用黑色小圆点标记出来
# 我们需要格点的经纬度坐标，并转换为地图投影下的坐标
lon2d, lat2d = np.meshgrid(sst_anomaly.longitude, sst_anomaly.latitude)
# 只选取significant为True的点
sig_lons = lon2d[significant]
sig_lats = lat2d[significant]

# 在地图上绘制散点，注意transform参数
ax.scatter(sig_lons, sig_lats, color='black', s=1, marker='.', alpha=0.6, transform=ccrs.PlateCarree())

ax.scatter用来打点，s=1控制点的大小，alpha=0.6让点有些透明度，不至于太突兀。这样，读者一眼就能看出哪些区域的暖或冷异常不是随机波动，而是具有统计意义的信号。

4.3 优化地图要素与输出

一张专业的图，细节决定成败。

# 优化海岸线和陆地
ax.add_feature(cfeature.COASTLINE.with_scale('50m'), linewidth=0.8, edgecolor='black')
ax.add_feature(cfeature.LAND, color='whitesmoke', edgecolor='none') # 陆地用更浅的颜色

# 添加国界线（根据绘图用途和规范决定是否添加）
# ax.add_feature(cfeature.BORDERS, linewidth=0.5, linestyle=':', edgecolor='gray')

# 使用CyclicColormap处理经度0/360度衔接问题（如果你的数据经度是0-360度）
# Cartopy的PlateCarree投影能自动处理，但填色时在0度经线附近可能有缝。
# 一个技巧是让数据经度从-180到180，或者使用cartopy的add_cyclic_point函数
from cartopy.util import add_cyclic_point
data_cyclic, lon_cyclic = add_cyclic_point(sst_anomaly.values, coord=sst_anomaly.longitude)
# 然后用data_cyclic和lon_cyclic去绘图

# 保存为矢量图，方便后期编辑和印刷
fig.savefig('global_sst_anomaly_2023_winter.pdf', dpi=300, bbox_inches='tight')
fig.savefig('global_sst_anomaly_2023_winter.svg', dpi=300, bbox_inches='tight')

使用add_cyclic_point可以避免在经度0度（即格林尼治子午线）附近出现一条难看的白色数据缝隙，这对于全球循环数据（如海温）的绘图是很好的实践。最后，保存为PDF或SVG矢量格式，无论放大多少倍都不会失真，在撰写论文时非常有用。

5. 常见问题排查与性能优化

即使按照步骤操作，你也可能会遇到一些“坑”。这里我总结几个最常见的问题和解决办法。

5.1 图形显示问题

问题1：地图是空白的，只有坐标轴。

• 检查：最可能的原因是ax.contourf或ax.pcolormesh中漏了transform=ccrs.PlateCarree()参数。Cartopy必须知道数据的坐标系。
• 检查：数据读取是否正确？用print(sst_anomaly.min(), sst_anomaly.max())看看数据是否有有效值（不是全0或全NaN）。
• 检查：地图范围set_extent是否设置正确，是否包含了你的数据范围？

问题2：颜色条显示不正常，全是同一个颜色。

• 检查：levels参数设置的范围是否远远大于或小于你数据的实际范围？比如数据在[-0.5, 0.5]之间，你却设置了levels=np.arange(-5,5,1)，那么所有数据都会落在同一个色阶区间内。调整levels或使用norm参数进行标准化。
• 检查：计算异常时，减法和平均的维度dim是否正确？确保是对time维度操作。

问题3：经度标签显示为0到360，但我想要-180到180。

• 处理：这是数据本身经度坐标的定义方式。你可以在绘图前转换数据坐标，也可以转换地图的显示方式。

  # 方法A：转换数据坐标（推荐，一劳永逸）
  ds = ds.assign_coords(longitude=(((ds.longitude + 180) % 360) - 180)).sortby('longitude')
  # 方法B：在绘图时，设置地图的中央经线为0，并设置extent为[-180,180]
  proj = ccrs.PlateCarree(central_longitude=0)
  ax.set_extent([-180, 180, -90, 90], crs=proj)
  

5.2 数据处理与性能

问题：数据文件太大，读取和计算速度慢。

• 策略1：延迟加载：xarray默认使用延迟加载，open_dataset并不会立刻把数据读入内存。只有在实际计算（如.mean()）或取值（如.values）时才会加载。利用这一点，先进行数据切片，再进行计算。

  # 不佳：先读入全部数据再切片
  # sst_all = ds.sst.load() # 这会立刻加载全部数据到内存
  # sst_winter = sst_all.where(ds.time.dt.month.isin([12,1,2]), drop=True)
  
  # 推荐：先切片，再计算（延迟计算链）
  sst_winter = ds.sst.sel(time=ds.time.dt.month.isin([12,1,2]))
  sst_winter_mean = sst_winter.mean(dim='time') # 此时才触发计算和加载
  

• 策略2：分块处理与并行计算：对于超大型数据集，可以使用xarray的chunk方法结合Dask库进行分块和并行计算。

  import dask.array as da
  # 以分块方式打开数据集
  ds_chunked = xr.open_dataset('big_data.nc', chunks={'time': 10, 'latitude': 100, 'longitude': 100})
  # 后续操作会自动并行化
  

• 策略3：降低空间分辨率：如果只是出图，并不需要原始数据的高分辨率（如0.25度）。可以在计算前对数据进行空间平均（coarsen或resample），比如降到1度分辨率，数据量会减少16倍，速度大幅提升。

5.3 图形美化与定制

问题：我想用其他投影，比如极射赤面投影看北极。

• 解决：Cartopy支持几十种地图投影。只需改变创建投影时的参数。

  # 北半球极射赤面投影
  proj_north = ccrs.NorthPolarStereo()
  ax_north = plt.axes(projection=proj_north)
  ax_north.set_extent([-180, 180, 60, 90], ccrs.PlateCarree()) # 设置显示北纬60度以上区域
  # 绘图时，transform参数仍然用PlateCarree，因为数据坐标没变
  cf = ax_north.contourf(lon, lat, data, transform=ccrs.PlateCarree(), ...)
  

问题：颜色条太宽或太窄，或者我想竖着放。

• 解决：调整plt.colorbar的参数。

  # 水平颜色条，放在图下方
  cbar = plt.colorbar(cf, ax=ax, orientation='horizontal', pad=0.05, shrink=0.8, aspect=30)
  # 垂直颜色条，放在图右侧
  cbar = plt.colorbar(cf, ax=ax, orientation='vertical', pad=0.05, shrink=0.6)
  # pad: 颜色条与地图的间距
  # shrink: 颜色条长度缩放因子
  # aspect: 颜色条的长宽比（针对水平条调节宽度）
  

画图本身是个不断调试和打磨的过程。我第一次画出来的图也是丑得没法看，经纬度标签叠在一起，颜色丑，海岸线粗糙。多试几次，多看看气象领域顶级期刊（如Journal of Climate）上的图是怎么做的，慢慢就能找到感觉。最关键的是理解每个参数背后的含义，这样遇到问题你才知道从哪里下手去调。希望这份详细的指南能帮你绕过我当年踩过的那些坑，顺利画出既科学又美观的专业气象图。