如何基于Matplotlib与Cartopy绘制多个圆形交集的轮廓边界

要实现多个测地圆的交集边界绘制，你需要借助shapely库的几何运算功能，结合Cartopy完成可视化。以下是修改后的完整代码：

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.geodesic import Geodesic  # 补全缺失的导入
from shapely.geometry import Polygon   # 补全缺失的导入
import numpy as np

fig = plt.figure(figsize=(10, 6))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.Mercator())

ax.set_extent([11, 38, 56, 70], crs=ccrs.PlateCarree())

ax.add_feature(cfeature.BORDERS, linestyle='-', linewidth=1)
ax.coastlines(linewidth=1)

# 添加圆形并收集所有多边形
radius = 120  # km
n_samples = 80  # 圆上的采样点数
lons = (21.64, 22.5, 23.82, 24.5, 25.44, 26.32, 26.9, 27.11, 27.38, 29.45, 29.8)
lats = (60.13, 61.81, 63.1, 60.56, 62.3, 64.77, 67.14, 60.9, 62.86, 63.84, 61.91)

# 收集所有圆形的Polygon对象
circle_polygons = []
for lon, lat in zip(lons, lats):
    circle_coords = Geodesic().circle(lon, lat, radius*1000., n_samples=n_samples)
    poly = Polygon(circle_coords)
    circle_polygons.append(poly)
    # 绘制原始圆形（保留你原来的可视化效果）
    feature = cfeature.ShapelyFeature([poly], ccrs.Geodetic(), fc='orange', alpha=0.4, ec="black", lw=0.5, linestyle="-")
    ax.add_feature(feature)

# 计算所有圆形的交集
if len(circle_polygons) > 0:
    intersection = circle_polygons[0]
    for poly in circle_polygons[1:]:
        intersection = intersection.intersection(poly)
        # 如果交集为空，提前终止计算
        if intersection.is_empty:
            break
    
    # 绘制交集的边界
    if not intersection.is_empty:
        # 处理交集可能是MultiPolygon的情况（多个不连通区域）
        if hasattr(intersection, 'geoms'):
            for geom in intersection.geoms:
                inter_feature = cfeature.ShapelyFeature([geom], ccrs.Geodetic(), ec='red', lw=2, fc='none', linestyle='-')
                ax.add_feature(inter_feature)
        else:
            inter_feature = cfeature.ShapelyFeature([intersection], ccrs.Geodetic(), ec='red', lw=2, fc='none', linestyle='-')
            ax.add_feature(inter_feature)

plt.show()

关键说明： #

• 补全依赖导入：原代码缺少Geodesic和Polygon的导入语句，这是代码运行的基础前提。
• 多边形收集：先将每个测地圆转换为shapely.Polygon对象并存储，为后续的几何运算做准备。
• 交集计算逻辑：从第一个多边形开始，依次与后续每个多边形求交集，逐步缩小范围得到所有圆的共同覆盖区域；如果中途交集为空，直接终止计算避免无效运算。
• 边界可视化：判断交集非空后，用红色粗线绘制交集边界（填充设为none仅显示轮廓），同时兼容交集为MultiPolygon的场景（即存在多个不连通的交集区域）。

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