在 NorthPolarStereo 投影中使用 Cartopy 绘制圆圈 [英] Drawing circles with Cartopy in NorthPolarStereo projection

问题描述

我想在 Cartopy 中以 NorthPolarStereo 投影绘制圆，以纬度、经度单位提供中心和半径.

Basemap 有类似的优秀问答

解决方案

圆心坐标必须是投影坐标.因此，需要进行坐标变换.

相关代码如下:

# 注意:lat = 72, lon = 100# proj = ccrs.NorthPolarStereo(central_longitude=lon)projx1, projy1 = proj.transform_point(lon, lat, ccrs.Geodetic()) #get proj coord of (lon,lat)ax.add_patch(mpatches.Circle(xy=[projx1, projy1], radius=r_ortho, color='red', \alpha=0.3，transform=proj，zorder=30))

输出图将是:

替代解决方案

由于使用的投影是conformal，天梭指示线图将始终是一个圆.这样，您需要的圆就可以用 Indicatrix 表示.这是您可以尝试的代码:

ax.tissot(rad_km=r_ortho/1000, lons=lon, lats=lat, n_samples=48, color='green', \alpha=0.3, zorder=31)

编辑 1

替换错误函数的代码:

导入pyprojdef compute_radius(val_degree):"""计算给定角度值(以度为单位)的表面距离(以米为单位)"""geod84 = pyproj.Geod(ellps='WGS84')lat0, lon0 = 0, 90_, _, dist_m = geod84.inv(lon0, lat0, lon0+val_degree, lat0)返回 dist_m计算半径(1)#获取:111319.49079327357(米)

I would like to draw circles in Cartopy in NorthPolarStereo projections, providing the center and radius in lat,lon units.

Similar and excellent questions and answers are available for Basemap here, and for Cartopy in Ortographic projection here. However, I would like to use the NorthPolarStereo in Cartopy. Trying the latter approach, just changing the projection makes the circle fixed in the North Pole, ignoring the coordinates you give for its center.

Any ideas on how to draw circles in Cartopy using NorthPolarStereo projection, providing its center and radius in lat,lon untis?

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 = 72 
lon = 100
r = 20

# Define the projection used to display the circle:
proj = ccrs.NorthPolarStereo(central_longitude=lon)


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)

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_patch(mpatches.Circle(xy=[lon, lat], radius=r_ortho, color='red', alpha=0.3, transform=proj, zorder=30))

plt.show()

Circle is fixed at North Pole, it won't move

解决方案

The coordinates of the circle's center must be the projection coordinates. So, a coordinate transformation is required.

Here is the relevant code:

# Note: lat = 72, lon = 100
# proj = ccrs.NorthPolarStereo(central_longitude=lon)

projx1, projy1 = proj.transform_point(lon, lat, ccrs.Geodetic()) #get proj coord of (lon,lat)
ax.add_patch(mpatches.Circle(xy=[projx1, projy1], radius=r_ortho, color='red', \
                            alpha=0.3, transform=proj, zorder=30))

The output plot will be:

Alternate Solution

Since the projection in use is conformal, a Tissot Indicatrix plot on it will always be a circle. So that, the circle you need can be represented with an Indicatrix. Here is the code you can try:

ax.tissot(rad_km=r_ortho/1000, lons=lon, lats=lat, n_samples=48, color='green', \
      alpha=0.3, zorder=31)

Edit 1

Code to replace the errorneous function:

import pyproj

def compute_radius(val_degree):
    """
    Compute surface distance in meters for a given angular value in degrees
    """
    geod84 = pyproj.Geod(ellps='WGS84')
    lat0, lon0 = 0, 90
    _, _, dist_m = geod84.inv(lon0, lat0,  lon0+val_degree, lat0)
    return dist_m

compute_radius(1)  # get: 111319.49079327357 (meters)

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