在Cartopy等距方位投影中使用contourf遇ValueError求助

问题概述 #

在Cartopy的Azimuthal Equidistant（等距方位）投影中使用contourf绘制全球信号传播时间时，出现以下错误：

line 185, in geos_multipolygon_from_polygons ValueError: Sequences of multi-polygons are not valid arguments

使用的Cartopy版本为0.21.1。

可复现代码 #

原始数据 #

# Longitudes of stations
longs = [-171.7827,  -171.7827,   179.1966,   179.1966,  -159.7733,  -159.7733,
  174.7043,   174.7043,   172.9229,   172.9229,   159.9475,   159.9475,
 -157.4457, -157.4457,   146.24998,  146.24998, -169.5292,  -169.5292,
  166.652,    166.652,   -155.5326,  -155.5326,  -158.0112,  -158.0112,
 -177.3698,  -177.3698,   144.8684,   166.7572,   166.7572,   117.239,
  117.239,    125.5791,   125.5791,   110.5354,   110.5354,   144.4382,
  144.4382,   138.20406,  138.20406, -176.6842,  -176.6842,   121.4971,
  121.4971,   126.62436,  126.62436,  -64.0489,   -64.0489,  -123.3046,
 -123.3046,  -110.7847,  -110.7847,   -90.2861,   -90.2861,  -106.4572,
 -106.4572,  -106.4572,  -147.8616,  -147.8616,  -147.8616,  -104.0359,
 -104.0359,   -95.83812,  -95.83812,  -70.7005,   -70.7005,    98.9443,
   98.9443,   -88.2763,   -88.2763,   -61.9787,  -61.9787,   -78.4508,
  -78.4508,  -175.385  ]

# Latitudes of stations
lats = [-13.9085,   -13.9085,    -8.5259,    -8.5259,   -21.2125,   -21.2125,
 -41.3087,   -41.3087,     1.3549,     1.3549,    -9.4387,    -9.4387,
   2.0448,     2.0448,   -20.08765,  -20.08765,   16.7329,    16.7329,
  19.2834,    19.2834,    19.7573,    19.7573,    21.42,      21.42,
  28.2156,    28.2156,    13.5893,   -77.8492,   -77.8492,   -32.9277,
 -32.9277,     7.0697,     7.0697,   -66.2792,   -66.2792,   -89.9289,
 -89.9289,    36.54567,   36.54567,   51.8823,    51.8823,    24.9735,
  24.9735,    37.47768,   37.47768, -64.7744,   -64.7744,    44.5855,
  44.5855,    32.3098,    32.3098,    -0.6742,    -0.6742,    34.94591,
  34.94591,   34.94591,   64.873599,  64.873599,  64.873599,  44.1212,
  44.1212,    29.96478,   29.96478,  -29.011,    -29.011,     18.8141,
  18.8141,    20.2263,    20.2263,   -38.0568,   -38.0568,     0.2376,
   0.2376,   -20.57    ]

# Time (h) signal detected after eruption
travel_time_h = [ 0.95296297,  0.95332528,  1.49046297,  1.4905475,   1.67046297, 
1.67026972, 2.3705475,   2.37046297,  2.60249194,  2.60240741,  2.7537963,   2.75360306,
  3.00943639,  3.00935186,  3.65610306,  3.65601852,  3.93165861,  3.93157408,
 16.13526972,  4.43074074,  4.61268519,  4.6130475,   4.6730475,   4.67296297,
  5.01026972,  5.01046297,  5.20768519,  5.96546297,  5.9655475,   6.49693639,
  6.49685186,  6.40324074,  6.40332528,  6.53740741,  6.53721417,  7.12074074,
  7.1205475,   7.34546297,  7.34499194,  7.26157408,  7.26221417,  7.64546297,
  7.6455475,   8.13407408,  8.13388083,  7.97693639,  7.97740741,  8.05082528,
  8.05101852,  8.00240741,  8.00221417,  8.65943639,  8.65907408,  8.41907408,
  8.41776972,  8.42722222,  8.94324074,  8.9430475,   8.94333333,  9.2555475,
  9.25601852,  8.99240741,  8.99249194,  9.26851852,  9.26749194,  9.16165861,
  9.16185186,  9.41990741,  9.41999194,  9.30851852,  9.31360306,  9.82324074,
  9.82332528,  0.        ]

原始插值与绘图代码 #

import matplotlib.pyplot as plt
import numpy as np
from cartopy import crs as ccrs
from scipy.interpolate import griddata

# Interpolate for contour
X, Y = np.meshgrid(longs, lats)
Z = griddata((longs, lats), travel_time_h, (X, Y), method='linear')

# Initialize figure
fig = plt.figure(figsize=(10, 8))
projLae = ccrs.AzimuthalEquidistant(central_longitude=-175.385, central_latitude=-20.57)
ax = plt.subplot(1, 1, 1, projection=projLae)

# Plot contour first as background
start_h, end_h, interval_h = 0.0, 10.0, 0.5
levels = np.arange(start=start_h, stop=end_h, step=interval_h)  # levels of contour
contour = ax.contourf(X, Y, Z, levels=levels, vmin=start_h, vmax=end_h, transform=ccrs.PlateCarree())

# Add colorbar for contour
cbar = fig.colorbar(contour, orientation='horizontal')   
cbar.ax.set_xlabel(f"Time [hr]")

# Plot station locations
ax.scatter(longs, lats, s=8, marker='*', color='red', transform=ccrs.PlateCarree())

# Plot map details
ax.coastlines()
ax.set_global()
plt.show()

错误原因分析 #

1. 网格生成问题：使用np.meshgrid(longs, lats)生成的是基于原始站点经纬度的不规则网格，这种网格在投影转换时会产生重叠或交叉的多边形，导致Cartopy的contourf无法正确处理。
2. 重复数据干扰：原始站点数据中存在大量重复的经纬度点，会影响插值结果的有效性，同时增加投影转换的复杂度。
3. 插值范围不足：原始插值仅覆盖站点所在的离散点区域，没有生成全球范围的规则网格，无法满足全球绘图需求。

修正方案 #

步骤1：生成规则的全球经纬度网格 #

替换原始插值部分，生成覆盖全球的规则网格，确保投影转换时的几何有效性。

步骤2：处理重复数据（可选） #

去除重复的站点数据，避免插值时的冗余计算。

步骤3：调整插值方法 #

使用cubic或nearest方法增强插值的平滑性，同时确保覆盖全球范围。

完整修正代码 #

import matplotlib.pyplot as plt
import numpy as np
from cartopy import crs as ccrs
from scipy.interpolate import griddata

# 去除重复的站点数据（可选，提升插值效率）
unique_indices = np.unique(np.array([longs, lats]).T, axis=0, return_index=True)[1]
unique_longs = np.array(longs)[unique_indices]
unique_lats = np.array(lats)[unique_indices]
unique_travel_time = np.array(travel_time_h)[unique_indices]

# 生成全球规则经纬度网格
lon_grid = np.linspace(-180, 180, 360)
lat_grid = np.linspace(-90, 90, 180)
X, Y = np.meshgrid(lon_grid, lat_grid)

# 插值生成全球传播时间数据
Z = griddata((unique_longs, unique_lats), unique_travel_time, (X, Y), method='cubic', fill_value=np.nan)
# 处理超出0-10小时的异常值
Z[Z > 10] = np.nan

# Initialize figure
fig = plt.figure(figsize=(10, 8))
projLae = ccrs.AzimuthalEquidistant(central_longitude=-175.385, central_latitude=-20.57)
ax = plt.subplot(1, 1, 1, projection=projLae)

# Plot contour first as background
start_h, end_h, interval_h = 0.0, 10.0, 0.5
levels = np.arange(start=start_h, stop=end_h, step=interval_h)  # levels of contour
contour = ax.contourf(X, Y, Z, levels=levels, vmin=start_h, vmax=end_h, 
                      transform=ccrs.PlateCarree(), cmap='viridis', extend='max')

# Add colorbar for contour
cbar = fig.colorbar(contour, orientation='horizontal', pad=0.05)   
cbar.ax.set_xlabel(f"Time [hr]")

# Plot station locations
ax.scatter(longs, lats, s=8, marker='*', color='red', transform=ccrs.PlateCarree(), zorder=5)

# Plot map details
ax.coastlines(linewidth=0.8)
ax.set_global()
plt.show()

关键修正点说明 #

• 生成规则全球网格：确保投影转换时的几何形状有效，避免无效多边形错误。
• 去除重复数据：减少冗余计算，提升插值稳定性。
• 处理异常值：将超出设定范围的16小时数据设为NaN，避免颜色映射异常。
• 添加zorder参数：确保站点标记显示在等高线之上。

内容的提问来源于stack exchange，提问作者just_another_profile