在正投影中用Cartopy绘制圆 [英] Drawing Circles with cartopy in orthographic projection

问题描述

我正在尝试使用Cartopy在给定的地理坐标下以一定的半径绘制圆.我想使用正投影法绘制，该投影法以圆心为中心.

I am trying to draw circles at a given geographical coordinate with a certain radius using cartopy. I want to draw using an orthographic projection, which is centred at the centre of the circle.

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

I use the following python code for testing:

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):

I get a correct result at the equator (set lat to 0 in example above):

但是当我朝杆子移动时，形状变形了(lat = 45):

But when I move towards a pole the shape is distorted (lat = 45):

在极点，我只看到圆的四分之一:

At the pole I only see one quarter of the circle:

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

I would expect to always see a perfect circle in orthographic projection, if the view is centred correctly. I also tried to use a different transform in the add_patch method, but then the circle completely vanishes!

推荐答案

您无法在PlateCarree坐标中定义圆，因为这是笛卡尔投影，并且在圆上绘制的圆不一定是球形的圆几何形状(除非您看到的圆是在(0，0)处).

You approach of defining the circle in PlateCarree coordinates is not going to work, because this is a cartesian projection and a circle drawn on it is not necessarily circular in spherical geometry (unless the circle is at (0, 0) as you saw).

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

Since you want the result to be circular in the Orthographic projection, you could draw the circle in native coordinates. This requires first defining your Orthographic projection centred on the centre of your circle, then computing what the radius of the circle (which you specify in degrees) would be in projection coordinates (distance from the centre of the projection). Doing it this way is convenient because it also gives you a neat way of determining the correct map extents. The example below computes the radius in orthographic coordinates by transforming a point 45 degrees north (or south if more convenient) away from the centre of the projection and gives the following:

完整代码如下:

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()

这篇关于在正投影中用Cartopy绘制圆的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！