python - 在正交投影中使用 cartopy 绘制圆
问题描述
我正在尝试使用 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):
如果视图正确居中,我希望始终在正交投影中看到一个完美的圆圈。我也尝试在add_patch
方法中使用不同的变换,但随后圆圈完全消失了!
解决方案
您在 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()
推荐阅读
- reactjs - 如何访问 React.js 中其他文件的数据?
- python - 如何使用 Python 和 mime 类型`text/plain` 编写一个字符文件?
- javascript - 切换容器后按钮未提交
- python - “系列”对象不是可调用的熊猫
- java - 在java中优化指数函数
- c - 使用指针的字符串长度不正确
- apache-kafka - 我无法在 kafka ksql 中创建表
- apache-flink - Flink:单独运行 sinks
- python - 如何在 Heroku 或类似设备上快速部署 Aiohttp Python 项目
- r - 在ggplot2中向图例添加简单的直方图