如何在Cartopy中正确绘制静止投影下的L2海表温度（SST）数据？

我尝试绘制L2海表温度（SST）数据，希望以静止投影在全球范围内展示，但输出结果不准确，部分陆地区域出现了不应存在的SST数值（表现为陆地区域被填充了SST颜色，与实际情况不符）。以下是我的代码：

import h5py
import sys
import numpy as np
import matplotlib.pyplot as plt    
import cartopy.crs as ccrs
import cartopy.feature as cfeature

# 从HDF5文件读取数据
fn = '/home/swadhin/project/insat/data/3RIMG_30MAR2018_0014_L2B_SST_V01R00.h5'

with h5py.File(fn) as f: 
    print(list(f.keys()))
    image = 'SST'
    img_arr = f[image][0,:,:]
    # 获取填充值用于数据掩码
    img_arr_fill = f[image].attrs['_FillValue'][0]
    # 从文件属性获取绘图范围
    left_lon = f.attrs['left_longitude'][0]
    right_lon = f.attrs['right_longitude'][0]
    lower_lat = f.attrs['lower_latitude'][0]
    upper_lat = f.attrs['upper_latitude'][0]
    sat_long = f.attrs['Nominal_Central_Point_Coordinates(degrees)_Latitude_Longitude'][1]
    sat_hght = f.attrs['Nominal_Altitude(km)'][0] * 1000.0 # 转换为米
print('Done reading HDF5 file')

# 用填充值掩码数据
img_arr_m = np.ma.masked_equal(img_arr, img_arr_fill)
print(img_arr_fill)
print(np.max(img_arr_m))
print(np.min(img_arr_m))

# 创建Geostationary投影绘图
map_proj = ccrs.Geostationary(central_longitude=sat_long,satellite_height=sat_hght)
ax = plt.axes(projection=map_proj)

ax.coastlines(color='black',linewidth = 0.5)
ax.gridlines(color='black', alpha=0.5, linestyle='--', linewidth=0.75, draw_labels=True)

map_extend_geos = ax.get_extent(crs=map_proj)
plt.imshow(img_arr_m, interpolation='none',origin='upper',extent=map_extend_geos, cmap = 'jet')
plt.colorbar()
plt.savefig('/home/swadhin/project/insat/data/l2_sst.png',format = 'png', dpi=1000)

我已准备好测试数据供调试使用。

问题根源 #

1. imshow的范围匹配错误：imshow默认将数据视为规则网格，直接使用ax.get_extent()获取的投影范围无法精准匹配SST数据的实际地理坐标，导致数据错位，陆地区域被错误填充。
2. 未处理陆地掩码：L2 SST数据中，陆地区域的数值通常不是_FillValue，而是有独立的掩码标识（比如质量标志、陆地掩码变量），仅掩码填充值不足以过滤陆地数据。

修改步骤及代码 #

1. 读取经纬度网格：从HDF文件中读取对应的经度和纬度数组，用于准确定位SST数据的地理坐标。
2. 用pcolormesh替代imshow：pcolormesh支持传入经纬度坐标和数据，能在Cartopy投影下更准确地渲染地理数据。
3. 添加陆地掩码处理：检查数据中的质量标志或陆地掩码变量，将陆地区域的SST值掩码；若无专用掩码，可显式绘制陆地特征覆盖异常区域。
4. 设置正确绘图范围：用数据的实际经纬度范围限制绘图区域，避免全球范围的无效渲染。

修改后的代码：

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

fn = '/home/swadhin/project/insat/data/3RIMG_30MAR2018_0014_L2B_SST_V01R00.h5'

with h5py.File(fn) as f: 
    print(list(f.keys()))
    # 读取SST数据
    sst_data = f['SST'][0, :, :]
    fill_value = f['SST'].attrs['_FillValue'][0]
    
    # 读取经纬度网格（根据HDF实际变量名调整）
    lon = f['Longitude'][:]
    lat = f['Latitude'][:]
    
    # 获取卫星参数
    sat_long = f.attrs['Nominal_Central_Point_Coordinates(degrees)_Latitude_Longitude'][1]
    sat_hght = f.attrs['Nominal_Altitude(km)'][0] * 1000.0
print('Done reading HDF5 file')

# 第一步：掩码填充值
sst_masked = np.ma.masked_equal(sst_data, fill_value)

# 第二步：掩码陆地区域（根据实际数据调整）
# 示例1：若有Land_Mask变量，1表示陆地
# land_mask = f['Land_Mask'][0, :, :]
# sst_masked = np.ma.masked_where(land_mask == 1, sst_masked)
# 示例2：若无专用掩码，用数值范围过滤异常值
# sst_masked = np.ma.masked_where((sst_masked > 50) | (sst_masked < -20), sst_masked)

# 设置投影
map_proj = ccrs.Geostationary(central_longitude=sat_long, satellite_height=sat_hght)
data_proj = ccrs.PlateCarree()  # 经纬度数据的坐标系

# 创建绘图对象
fig, ax = plt.subplots(figsize=(12, 10), subplot_kw={'projection': map_proj})

# 添加地图特征
ax.coastlines(color='black', linewidth=0.5)
ax.add_feature(cfeature.LAND, facecolor='lightgray')  # 显式绘制陆地，覆盖异常区域
ax.gridlines(color='black', alpha=0.5, linestyle='--', linewidth=0.75, draw_labels=True)

# 绘制SST数据，指定数据坐标系
mesh = ax.pcolormesh(lon, lat, sst_masked, transform=data_proj, cmap='jet', shading='auto')

# 添加颜色条
plt.colorbar(mesh, ax=ax, orientation='horizontal', pad=0.05, label='SST (°C)')

# 设置绘图范围为数据实际经纬度范围
ax.set_extent([lon.min(), lon.max(), lat.min(), lat.max()], crs=data_proj)

# 保存图像
plt.savefig('/home/swadhin/project/insat/data/l2_sst_fixed.png', format='png', dpi=300)
plt.close()

关键说明 #

• 经纬度网格：确保HDF文件中存在Longitude和Latitude变量，若变量名不同需对应修改。
• 陆地掩码：优先使用数据自带的陆地掩码变量；若无，可通过SST数值范围过滤，或显式添加cfeature.LAND覆盖陆地区域。
• 投影匹配：pcolormesh的transform参数必须设置为数据的坐标系（通常是PlateCarree），让Cartopy自动转换到Geostationary投影。

内容的提问来源于stack exchange，提问作者The Emerging Star