如何用Matplotlib在不同投影下绘制Tissot's indicatrix？

你的问题核心在于：当前代码是在**投影平面（PlateCarree）**上绘制固定半径的圆，而非先在球面上定义小圆再做投影转换，因此无法体现投影带来的拉伸/压缩变形。Tissot指示线本质是球面上的等半径小圆，必须先在地理坐标系（球面）定义，再通过投影转换到目标坐标系，才能呈现出正确的变形效果。

修改后的代码 #

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import math

def create_tissot_indicatrix(center_lon, center_lat, radius_deg, num_points=100):
    # 生成球面上指定半径的小圆经纬度点
    geod = ccrs.Geodetic()
    lon_points, lat_points = [], []
    # 遍历角度生成圆周上的点
    for angle in range(num_points):
        theta = math.radians(angle * 360 / num_points)
        # 利用地理坐标系正算方法，计算从中心点出发的目标点
        # 1度纬度对应的地面距离约为111319.9米
        lon, lat, _ = geod.fwd(center_lon, center_lat, math.degrees(theta), radius_deg * 111319.9)
        lon_points.append(lon)
        lat_points.append(lat)
    # 闭合路径确保图形完整
    lon_points.append(lon_points[0])
    lat_points.append(lat_points[0])
    return lon_points, lat_points

def main():
    # 可切换PlateCarree或Mercator投影测试
    proj = ccrs.PlateCarree()
    # proj = ccrs.Mercator()
    
    lonW = -90
    lonE = 90
    latS = -90
    latN = 90
    res = '110m'

    fig = plt.figure(figsize=(11, 8.5))
    ax = plt.subplot(1, 1, 1, projection=proj)
    ax.gridlines(draw_labels=True, linewidth=2, color='gray', alpha=0.5, linestyle='--')
    ax.set_extent([lonW, lonE, latS, latN], crs=proj)
    ax.coastlines(resolution=res, color='black')

    ax.add_feature(cfeature.LAND)
    ax.add_feature(cfeature.OCEAN)
    ax.add_feature(cfeature.BORDERS, linestyle='-', linewidth=0.5)
    ax.add_feature(cfeature.LAKES)

    # 添加单个坐标点
    ax.scatter(24.3, 61.83, color='red', marker='o', s=6, transform=ccrs.Geodetic(), label='Hyytiälä')

    # 绘制不同纬度的Tissot指示线
    # 赤道
    lon_points, lat_points = create_tissot_indicatrix(0, 0, 10)
    ax.fill(lon_points, lat_points, color='red', alpha=0.3, transform=ccrs.Geodetic(), zorder=30)
    
    # 北纬30度
    lon_points, lat_points = create_tissot_indicatrix(0, 30, 10)
    ax.fill(lon_points, lat_points, color='blue', alpha=0.3, transform=ccrs.Geodetic(), zorder=30)
    
    # 北纬60度
    lon_points, lat_points = create_tissot_indicatrix(0, 60, 10)
    ax.fill(lon_points, lat_points, color='green', alpha=0.3, transform=ccrs.Geodetic(), zorder=30)
    
    # 北极点
    lon_points, lat_points = create_tissot_indicatrix(0, 90, 10)
    ax.fill(lon_points, lat_points, color='purple', alpha=0.3, transform=ccrs.Geodetic(), zorder=30)

    plt.legend()
    plt.show()

if __name__ == "__main__":
    main()

关键说明 #

1. 地理坐标系转换：使用ccrs.Geodetic()作为transform参数，告诉Cartopy所有点都是球面地理坐标，会自动转换到当前的投影坐标系（PlateCarree或Mercator）。
2. 球面小圆生成：通过geod.fwd()方法计算球面上的圆周点，确保生成的是球面上半径固定的小圆，而非投影平面上的圆。
3. 投影切换测试：只需修改proj变量为ccrs.Mercator()，即可看到Mercator投影下靠近极点的指示线被明显拉伸的效果，完全符合Tissot指示线的特性。

内容的提问来源于stack exchange，提问作者Louschmuh