如何用cartopy绘制跨越极点的多边形

问题描述 投票:0回答:1

我试图在地图上的任意位置绘制多边形,包括跨越极点和日期变更线的位置。从概念上讲,考虑绘制用于轨道测量的仪器足迹,其中足迹的角以纬度/经度(或笛卡尔)表示。我已经能够绘制中纬度多边形,但跨越极点的形状却显示为空白。

这里有一些代码来显示我得到的结果:

from matplotlib.patches import Polygon
import matplotlib.pyplot as plt
import cartopy.crs as ccrs

# https://stackoverflow.com/questions/73689586/
geodetic = ccrs.Geodetic()

# Define a polygon over the pole (long, lat)
points_north = [ (0,   75),
                 (270, 75),
                 (180, 75),
                 (90,  75), ]

# Mid-latitude polygon (long, lat)
points_mid = [ (10,  70),
               (30,  20),
               (60,  20),
               (80,  70), ]

# Create a PlateCarree projection
proj0 = ccrs.PlateCarree()
proj0._threshold /= 5.                 # https://stackoverflow.com/questions/59020032/
ax = plt.axes( projection = proj0 )

# Add the polygons to the plot
ax.add_patch( Polygon( points_north, alpha = 0.2, transform = geodetic ) )
ax.add_patch( Polygon( points_mid,   alpha = 0.2, transform = geodetic ) )

# Some window dressing
longlocs = list( range( -180, 181, 30 ) )
latlocs = [ -75, -60, -30, 0, 30, 60, 75 ]
ax.gridlines( xlocs = longlocs, ylocs = latlocs,
              draw_labels = True,  crs = proj0 )
ax.set_global()
ax.set_xlabel('Longitude')
ax.set_ylabel('Latitude')

# Show the plot
plt.show()

这会产生以下地图:

多边形投影

中纬度多边形按预期绘制,但北极正方形没有。我希望该多边形能够以扇形边缘涂抹在整个北部,类似于旋转极盒示例中所示的内容。

我已经尝试过 RotatedPole,但我真的还没有接近一个可以让我绘制任意足迹的通用解决方案。

python matplotlib cartopy
1个回答
0
投票

PlateCarree
投影有以下不规则之处:-

  1. 北/南极退化为线,1点-->无穷大点
  2. 东西经线180度分成2条不同的线(实际上它们是同一条线)

任何其部分包含/穿过退化的几何体也可能继承此类不规则性。

根据这些限制,无法绘制覆盖 N 极和跨越日期变更线的多边形。

为了避免“退化问题”,您需要将多边形切割成几块,并且每块都必须小心地设置在退化区域之外。

这里是与您的

points_north
等效的4组点,每个点都有90度的经度范围,由前2个点和N极的第3个点定义,以完成多边形边缘。

tri1 = [(-180, 75), (-90, 75), (-90, 90)]
tri2 = [(-90, 75), (0, 75), (0, 90)]
tri3 = [(0, 75), (90, 75), (90, 90)]
tri4 = [(90, 75), (180, 75), (180, 90)]

修改后的代码:-

from matplotlib.patches import Polygon
import matplotlib.pyplot as plt
import cartopy.crs as ccrs

geodetic = ccrs.Geodetic()

#points_north = [ (0,  75),
#                (270, 75),
#                (180, 75),
#                (90,  75)]

# New 4 sets of points that can create (triangle) polygons
tri1 = [(-180, 75), (-90, 75), (-90, 90)]
tri2 = [(-90, 75), (0, 75), (0, 90)]
tri3 = [(0, 75), (90, 75), (90, 90)]
tri4 = [(90, 75), (180, 75), (180, 90)]

# Mid-latitude polygon (long, lat)
points_mid = [ (10,  70),
               (30,  20),
               (60,  20),
               (80,  70)]

proj0 = ccrs.PlateCarree()
#proj0 = ccrs.Orthographic(central_longitude=-30, central_latitude=60)
proj0._threshold /= 10
ax = plt.axes( projection = proj0 )

# Add the polygons to the plot
ax.add_patch( Polygon( tri1, alpha = 0.7, lw=1, transform = geodetic ) )
ax.add_patch( Polygon( tri2, alpha = 0.7, lw=1, transform = geodetic ) )
ax.add_patch( Polygon( tri3, alpha = 0.7, lw=1, transform = geodetic ) )
ax.add_patch( Polygon( tri4, alpha = 0.7, lw=1, transform = geodetic ) )

ax.add_patch( Polygon( points_mid, lw=2, ec='red', fc='gray', alpha=0.35, transform=geodetic ) )

ax.gridlines( draw_labels = True,  crs = ccrs.PlateCarree() )
ax.set_global()
plt.show()

输出图: fig1

更合适的投影是

Orthographic
投影。 你可以尝试使用这个:-

proj0 = ccrs.Orthographic(central_longitude=-30, central_latitude=60)

代替 PlateCarree 投影,并得到此图。

fig2

© www.soinside.com 2019 - 2024. All rights reserved.