Cartopy绘制大圈出现分裂问题求助

问题原因 #

当cgeod.Geodesic().circle()生成的圆形多边形跨越180°经线（或投影的中央经线边界）时，Shapely多边形的点列顺序会导致Matplotlib的填充逻辑出错，表现为图形沿纬线分裂并填充了错误的区域（实际填充了圆形外部而非内部）。

解决方案 #

方法1：调整投影中央经线，避免多边形跨边界 #

将投影的中央经线设置为丹佛的经度，让生成的圆形完全处于投影的可视范围内，避免跨边界导致的填充异常：

cm = -104.985  # 使用丹佛的经度作为投影中央经线
proj = crs.PlateCarree(central_longitude=cm)

方法2：修复多边形并选择正确的填充区域 #

若必须使用中央经线为0的投影，可通过修复多边形环方向、判断面积来选择需要填充的圆形内部区域：

# 生成圆形点后处理多边形
circle_points = cgeod.Geodesic().circle(lon=site_lon, lat=site_lat, radius=dFOV, n_samples=1000, endpoint=False)
# 修复无效的多边形环
geom = shapely.geometry.Polygon(circle_points).buffer(0)
# 筛选出面积较小的区域（即圆形内部，大的是地球剩余区域）
if geom.area > 1e10:
    full_earth = shapely.geometry.box(-180, -90, 180, 90)
    geom = full_earth.difference(geom)

完整修正代码（方法1示例） #

import numpy as np
import cartopy.geodesic as cgeod
import cartopy.crs as crs
import cartopy.feature as cfeature
import shapely
import matplotlib.pyplot as plt

# SET CONSTANTS
RE = 6371008.8 # meters, WGS84 Sphere
h = 35786000.0 # meters, GEO

def arc_dist_to_FOV(h, ah):
    a = (np.pi/180)*(90 + ah) # radians
    b = np.arcsin(RE * np.sin(a)/(RE + h)) # radians
    g = np.pi - a - b # radians
    d = (180/np.pi)*g*RE
    return d

# 用丹佛经度作为中央经线，避免多边形跨投影边界
cm = -104.985
proj = crs.PlateCarree(central_longitude=cm)
fig = plt.figure(figsize=(10,8))
ax = fig.add_subplot(1,1,1, projection=proj)

ax.set_global()

ax.add_feature(cfeature.COASTLINE, edgecolor='black')
ax.add_feature(cfeature.BORDERS, edgecolor='black')
ax.add_feature(cfeature.STATES, edgecolor='black')

ax.stock_img()
ax.gridlines(draw_labels=True, crs=proj)
lat, lon = 39.7392, -104.985

plt.scatter(x=lon, y=lat, color='blue', s=10, transform=proj)

# COMPUTE FOV
ah = 20 # degrees above horizon
dFOV = arc_dist_to_FOV(h, ah)
site_lat = lat
site_lon = lon
    
# ADD SHAPES TO MAP
circle_points = cgeod.Geodesic().circle(lon=site_lon, lat=site_lat, radius=dFOV, n_samples=1000, endpoint=False)
geom = shapely.geometry.Polygon(circle_points)
ax.add_geometries((geom,), crs=proj, alpha=0.5, facecolor='red', edgecolor='black', linewidth=1)

plt.show()

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