在正交投影中用 cartopy 绘制圆圈

Question

我正在尝试使用 cartopy 在给定的地理坐标处以一定的半径绘制圆。我想使用以圆心为中心的正交投影进行绘制。

我使用以下python代码进行测试：

import numpy as np
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches

# example: draw circle with 45 degree radius around the North pole
lon = 0
lat = 90
r = 45

# find map ranges (with 5 degree margin)
minLon = lon - r - 5
maxLon = lon + r + 5
minLat = lat - r - 5
maxLat = lat + r + 5

# define image properties
width = 800
height = 800
dpi = 96
resolution = '50m'

# create figure
fig = plt.figure(figsize=(width / dpi, height / dpi), dpi=dpi)
ax = fig.add_subplot(1, 1, 1, projection=ccrs.Orthographic(central_longitude=lon, central_latitude=lat))
ax.set_extent([minLon, maxLon, minLat, maxLat])
ax.imshow(np.tile(np.array([[cfeature.COLORS['water'] * 255]], dtype=np.uint8), [2, 2, 1]), origin='upper', transform=ccrs.PlateCarree(), extent=[-180, 180, -180, 180])
ax.add_feature(cfeature.NaturalEarthFeature('physical', 'land', resolution, edgecolor='black', facecolor=cfeature.COLORS['land']))
ax.add_feature(cfeature.NaturalEarthFeature('cultural', 'admin_0_countries', resolution, edgecolor='black', facecolor='none'))
ax.add_feature(cfeature.NaturalEarthFeature('physical', 'lakes', resolution, edgecolor='none', facecolor=cfeature.COLORS['water']), alpha=0.5)
ax.add_feature(cfeature.NaturalEarthFeature('physical', 'rivers_lake_centerlines', resolution, edgecolor=cfeature.COLORS['water'], facecolor='none'))
ax.add_feature(cfeature.NaturalEarthFeature('cultural', 'admin_1_states_provinces_lines', resolution, edgecolor='gray', facecolor='none'))
ax.add_patch(mpatches.Circle(xy=[lon, lat], radius=r, color='red', alpha=0.3, transform=ccrs.PlateCarree(), zorder=30))
fig.tight_layout()
plt.savefig('CircleTest.png', dpi=dpi)

plt.show()

我在赤道处得到了正确的结果（在上面的示例中将 lat 设置为 0）：

但是当我向一根杆子移动时，形状会变形 (lat = 45)：

在极点我只看到了圆的四分之一：

如果视图居中正确，我希望在正交投影中总是看到一个完美的圆。我还尝试在 add_patch 方法中使用不同的变换，但是圆圈完全消失了！

Answer 1

你在 PlateCarree 坐标中定义圆的方法是行不通的，因为这是一个笛卡尔投影，在它上面绘制的圆在球面几何中不一定是圆形（除非圆在 (0, 0)如你所见）。

由于您希望结果在正交投影中是圆形的，因此您可以在本地坐标中绘制圆形。这需要首先定义以圆心为中心的正交投影，然后计算投影坐标（距投影中心的距离）中圆的半径（以度为单位指定）。这样做很方便，因为它还为您提供了一种确定正确地图范围的简洁方法。下面的示例通过将一个点从投影中心向北（或向南，如果更方便的话）转换 45 度来计算正交坐标中的半径，并给出以下内容：

完整代码如下：

import numpy as np
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches

# example: draw circle with 45 degree radius around the North pole
lat = 51.4198101
lon = -0.950854653584
r = 45

# Define the projection used to display the circle:
proj = ccrs.Orthographic(central_longitude=lon, central_latitude=lat)


def compute_radius(ortho, radius_degrees):
    phi1 = lat + radius_degrees if lat <= 0 else lat - radius_degrees
    _, y1 = ortho.transform_point(lon, phi1, ccrs.PlateCarree())
    return abs(y1)

# Compute the required radius in projection native coordinates:
r_ortho = compute_radius(proj, r)

# We can now compute the correct plot extents to have padding in degrees:
pad_radius = compute_radius(proj, r + 5)

# define image properties
width = 800
height = 800
dpi = 96
resolution = '50m'

# create figure
fig = plt.figure(figsize=(width / dpi, height / dpi), dpi=dpi)
ax = fig.add_subplot(1, 1, 1, projection=proj)
# Deliberately avoiding set_extent because it has some odd behaviour that causes
# errors for this case. However, since we already know our extents in native
# coordinates we can just use the lower-level set_xlim/set_ylim safely.
ax.set_xlim([-pad_radius, pad_radius])
ax.set_ylim([-pad_radius, pad_radius])
ax.imshow(np.tile(np.array([[cfeature.COLORS['water'] * 255]], dtype=np.uint8), [2, 2, 1]), origin='upper', transform=ccrs.PlateCarree(), extent=[-180, 180, -180, 180])
ax.add_feature(cfeature.NaturalEarthFeature('physical', 'land', resolution, edgecolor='black', facecolor=cfeature.COLORS['land']))
ax.add_feature(cfeature.NaturalEarthFeature('cultural', 'admin_0_countries', resolution, edgecolor='black', facecolor='none'))
ax.add_feature(cfeature.NaturalEarthFeature('physical', 'lakes', resolution, edgecolor='none', facecolor=cfeature.COLORS['water']), alpha=0.5)
ax.add_feature(cfeature.NaturalEarthFeature('physical', 'rivers_lake_centerlines', resolution, edgecolor=cfeature.COLORS['water'], facecolor='none'))
ax.add_feature(cfeature.NaturalEarthFeature('cultural', 'admin_1_states_provinces_lines', resolution, edgecolor='gray', facecolor='none'))
ax.add_patch(mpatches.Circle(xy=[lon, lat], radius=r_ortho, color='red', alpha=0.3, transform=proj, zorder=30))
fig.tight_layout()
plt.savefig('CircleTest.png', dpi=dpi)
plt.show()

Answer 2

这可能有点晚了，但是 Cartopy 中有一个方便的功能。

我们可以使用 Cartopy 的 .circle 函数 (documentation) 从测地坐标系中的特定（经度和纬度）生成具有指定半径的点环，然后用这些点绘制多边形使用 Shapely.

这看起来像下面这样

circle_points = cartopy.geodesic.Geodesic().circle(lon=lon, lat=lat, radius=radius_in_meters, n_samples=n_points, endpoint=False)
geom = shapely.geometry.Polygon(circle_points)
ax.add_geometries((geom,), crs=cartopy.crs.PlateCarree(), facecolor='red', edgecolor='none', linewidth=0)

将 crs 指定为 PlateCarree 并不重要，只是避免了 Shapely 的警告。你会保持你想要的投影。然而，如果你直接用圆心在极点上绘图，你可能仍然有问题，可能需要做一些花哨的转换（最近没有测试过，但回想几个月前它有点不稳定） .

您也可以使用 Cartopy 使用的 pyproj 库手动计算这些点，特别是 Geod class。选择一个具有半径的点并通过 azmoths 循环，无论你希望你的圆与 .inv 或 .fwd 函数的精细程度类似于中的建议。我不推荐这种方法，但很久以前用它来完成同样的事情。

在正交投影中用 cartopy 绘制圆圈

Drawing Circles with cartopy in orthographic projection

python

orthographic

cartopy