如何使用Cartopy和Matplotlib绘制Tissot？

Question

如何使用Cartopy和Matplotlib绘制Tissot？

matplotlibcartopy

5

我刚刚从Basemap转换到cartopy，用于绘制星空图，我更喜欢它。

(主要原因是在某些电脑上Basemap出现故障，我无法修复。)

我唯一遇到的问题是如何得到一个Tissot圆（用于显示望远镜的视野锥形）。

这是一些绘制随机星星的示例代码（我在真实场景中使用目录）：

import matplotlib.pyplot as plt
from cartopy import crs
import numpy as np

# create some random stars:

n_stars = 100
azimuth = np.random.uniform(0, 360, n_stars)
altitude = np.random.uniform(75, 90, n_stars)
brightness = np.random.normal(8, 2, n_stars)

fig = plt.figure()
ax = fig.add_subplot(1,1,1, projection=crs.NorthPolarStereo())
ax.background_patch.set_facecolor('black')

ax.set_extent([-180, 180, 75, 90], crs.PlateCarree())

plot = ax.scatter(
    azimuth,
    altitude,
    c=brightness,
    s=0.5*(-brightness + brightness.max())**2,
    transform=crs.PlateCarree(),
    cmap='gray_r',
)

plt.show()

我该如何在该图像中添加一个特定半径的Tissot圆（以度为单位）？ https://en.wikipedia.org/wiki/Tissot%27s_indicatrix

- MaxNoe

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- pelson · Accepted Answer

我一直打算回去并添加GeographicLib中提供正反解大圆计算的两个函数，这只需要在给定纬度/经度/半径的情况下通过适当的方位角进行采样就可以计算出一个大地圆。可惜，我还没有做到这一点，但是在pyproj中有一个相当原始（但有效）的包装器来实现此功能。

要实现一个Tissot指数，代码可能会如下所示：

import matplotlib.pyplot as plt

import cartopy.crs as ccrs
import numpy as np

from pyproj import Geod
import shapely.geometry as sgeom


def circle(geod, lon, lat, radius, n_samples=360):
    """
    Return the coordinates of a geodetic circle of a given
    radius about a lon/lat point.

    Radius is in meters in the geodetic's coordinate system.

    """
    lons, lats, back_azim = geod.fwd(np.repeat(lon, n_samples),
                                     np.repeat(lat, n_samples),
                                     np.linspace(360, 0, n_samples),
                                     np.repeat(radius, n_samples),
                                     radians=False,
                                     )
    return lons, lats


def main():
    ax = plt.axes(projection=ccrs.Robinson())
    ax.coastlines()

    geod = Geod(ellps='WGS84')

    radius_km = 500
    n_samples = 80

    geoms = []
    for lat in np.linspace(-80, 80, 10):
        for lon in np.linspace(-180, 180, 7, endpoint=False):
            lons, lats = circle(geod, lon, lat, radius_km * 1e3, n_samples)
            geoms.append(sgeom.Polygon(zip(lons, lats)))

    ax.add_geometries(geoms, ccrs.Geodetic(), facecolor='blue', alpha=0.7)

    plt.show()


if __name__ == '__main__':
    main()