使用Matplotlib+contourf绘制全球数据时遇拓扑错误求助

尝试使用静止轨道投影（Geostationary projection）绘制全球范围数据时，出现TopologicalError，推测与边界区域有关。相关代码、报错信息如下：

运行代码 #

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 = '..../3DIMG_03MAY2019_0230_L1B_STD_V01R00.h5'

with h5py.File(fn) as f:
    # 获取影像数据
    image = 'IMG_TIR1'
    image2 = 'IMG_TIR2'
    lon = f['Longitude'][:]*0.01
    lat = f['Latitude'][:]*0.01
    img_arr = f[image][0, :, :]
    img_arr2 = f[image2][0, :, :]
    
    # 获取填充值用于数据掩码
    img_arr_fill = f[image].attrs['_FillValue'][0]
    img_arr_fill2 = f[image2].attrs['_FillValue'][0]

    sat_long = f.attrs['Nominal_Central_Point_Coordinates(degrees)_Latitude_Longitude'][1]
    sat_hght = f.attrs['Nominal_Altitude(km)'][0] * 1000.0  # 转换为米
    bt_lut_tir1 = np.array(f[image+str('_TEMP')])
    bt_lut_tir2 = np.array(f[image2+str('_TEMP')])
print('HDF5文件读取完成')

# 用np.ma.masked_equal掩码掉角落的填充值数据
img_arr_m = np.ma.masked_equal(img_arr, img_arr_fill)
img_arr_m2 = np.ma.masked_equal(img_arr2, img_arr_fill2)
lon_m = np.ma.masked_equal(lon, 327.67)
lat_m = np.ma.masked_equal(lat, 327.67)

# 通过查找表计算亮温
def count2bt_lut_tir1(count):
    return bt_lut_tir1[count]
def count2bt_lut_tir2(count):
    return bt_lut_tir2[count]
bt_tir1_lut = count2bt_lut_tir1(img_arr_m)
bt_tir2_lut = count2bt_lut_tir2(img_arr_m2)

######### 绘图部分 
# 使用cartopy和matplotlib创建静止轨道投影图  

map_proj = ccrs.Geostationary(central_longitude=sat_long,satellite_height=sat_hght)

ax = plt.axes(projection=map_proj)
ax.set_global()
ax.coastlines(color='black',linewidth = 0.5)
ax.add_feature(cfeature.BORDERS, edgecolor='white', linewidth=0.25)
ax.add_feature(cfeature.STATES,edgecolor = 'red',linewidth = 0.5)
ax.gridlines(color='black', alpha=0.5, linestyle='--', linewidth=0.75, draw_labels=True)
plt.title('Brightness Temperature(K) in TIR1 (10.8 micron)')
cb = ax.contourf(lon_m,lat_m,bt_tir1_lut,cmap = 'jet',transform = ccrs.PlateCarree())
plt.colorbar(cb)
plt.show()

报错信息 #

---------------------
PredicateError                            Traceback (most recent call last)
File ~/anaconda3/envs/rttov/lib/python3.9/site-packages/shapely/predicates.py:15, in BinaryPredicate.__call__(self, this, other, *args)
     14 try:
---> 15     return self.fn(this._geom, other._geom, *args)
     16 except PredicateError as err:
     17     # Dig deeper into causes of errors.

File ~/anaconda3/envs/rttov/lib/python3.9/site-packages/shapely/geos.py:584, in errcheck_predicate(result, func, argtuple)
    583 if result == 2:
--> 584     raise PredicateError("Failed to evaluate %s" % repr(func))
    585 return result

PredicateError: Failed to evaluate <_FuncPtr object at 0x7fb67e380f40>

During handling of the above exception, another exception occurred:

TopologicalError                          Traceback (most recent call last)
File ~/anaconda3/envs/rttov/lib/python3.9/site-packages/IPython/core/formatters.py:339, in BaseFormatter.__call__(self, obj)
    337     pass
    338 else:
--&gt; 339     return printer(obj)
    340 # Finally look for special method names
    341 method = get_real_method(obj, self.print_method)
...
     37             "Likely cause is invalidity of the geometry %s" % (
     38                 self.fn.__name__, repr(geom)))
     39 raise err

TopologicalError: The operation 'GEOSContains_r' could not be performed. Likely cause is invalidity of the geometry <shapely.geometry.polygon.Polygon object at 0x7fb5c0ae45b0>

这个TopologicalError通常是Cartopy加载的地理要素（如边界、州界）存在无效几何图形，或数据投影转换时触发边界校验问题，可尝试以下修复方式：

• 临时禁用可疑地理要素：
  报错指向无效几何图形，大概率是cfeature.STATES或cfeature.BORDERS的内置数据有问题。先注释掉这两行代码，测试是否能正常运行：

  # ax.add_feature(cfeature.BORDERS, edgecolor='white', linewidth=0.25)
  # ax.add_feature(cfeature.STATES,edgecolor = 'red',linewidth = 0.5)
  

  若运行正常，后续可替换为第三方可靠shapefile数据源，而非默认cfeature。

• 转换数据投影后再绘制：
  静止轨道卫星数据本身基于星下点投影，直接用经纬度绘制易在边缘触发投影冲突。可先将经纬度转换为静止轨道投影坐标：

  # 将经纬度转换为Geostationary投影坐标
  coords = map_proj.transform_points(ccrs.PlateCarree(), lon_m, lat_m)
  x = coords[..., 0]
  y = coords[..., 1]
  # 绘制时无需指定transform参数
  cb = ax.contourf(x, y, bt_tir1_lut, cmap='jet')
  

• 修复无效几何图形：
  若需保留地理要素，可借助shapely工具修复无效几何：

  from shapely.geometry import Polygon
  from shapely.validation import make_valid
  
  # 自定义修复地理要素的函数
  def get_valid_feature(feature):
      valid_geoms = []
      for geom in feature.geometries():
          if not geom.is_valid:
              valid_geom = make_valid(geom)
              if isinstance(valid_geom, Polygon):
                  valid_geoms.append(valid_geom)
              else:
                  valid_geoms.extend([g for g in valid_geom.geoms if isinstance(g, Polygon)])
          else:
              valid_geoms.append(geom)
      return cfeature.ShapelyFeature(valid_geoms, feature.crs)
  
  # 使用修复后的要素
  borders = get_valid_feature(cfeature.BORDERS)
  ax.add_feature(borders, edgecolor='white', linewidth=0.25)
  
  states = get_valid_feature(cfeature.STATES)
  ax.add_feature(states, edgecolor='red', linewidth=0.5)
  

• 缩小绘图范围：
  静止轨道卫星观测范围有限，ax.set_global()强制绘制全球会让边缘无效数据触发几何校验。可根据卫星星下点设置合理范围：

  # 设置显示±60度经度范围
  ax.set_extent([sat_long-60, sat_long+60, -60, 60], crs=ccrs.PlateCarree())
  

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