使用Cartopy绘制等距方位投影大圆路径与指定坐标错位问题

问题描述 #

使用Python的Cartopy和Matplotlib库，在等距方位投影（Azimuthal Equidistant projection）上绘制开普敦与霍巴特两点间的大圆路径时，代表两点的标记位置正确，但生成的大圆路径出现偏移，无法准确连接两点。

原代码如下：

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

def great_circle_points(lat1, lon1, lat2, lon2, num_points):
    lat1, lon1, lat2, lon2 = np.radians([lat1, lon1, lat2, lon2])
    dlat = lat2 - lat1
    dlon = lon2 - lon1
    fractions = np.linspace(0, 1, num_points)
    a = np.sin(dlat/2)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2)**2
    delta_sigma = 2 * np.arctan2(np.sqrt(a), np.sqrt(1-a))

    lat = np.arctan2(
        np.sin(lat1) * np.cos(delta_sigma * fractions) + np.sin(lat2) * np.sin(delta_sigma * fractions),
        np.sqrt((np.cos(lat1) + np.cos(lat2) * np.cos(delta_sigma * fractions))**2 + (np.cos(lat2) * np.sin(delta_sigma * fractions))**2)
    )
    lon = lon1 + np.arctan2(
        np.sin(delta_sigma * fractions) * np.sin(lon2 - lon1),
        np.cos(lat1) * np.cos(delta_sigma * fractions) - np.sin(lat1) * np.sin(delta_sigma * fractions) * np.cos(lon2 - lon1)
    )

    lat = np.degrees(lat)
    lon = np.degrees(lon)

    return lat, lon

lat1, lon1 = -33.9715, 18.6021  # Cape Town
lat2, lon2 = -42.8364, 147.5075  # Hobart
num_points = 1001

latitudes, longitudes = great_circle_points(lat1, lon1, lat2, lon2, num_points)

plt.figure(figsize=[10, 10])
#ccrs.Orthographic(), ccrs.Stereographic(), ccrs.LambertAzimuthalEqualArea(), or ccrs.Gnomonic().
ax = plt.axes(projection=ccrs.AzimuthalEquidistant(central_latitude=90, central_longitude=0))
ax.set_global()
ax.coastlines()

# Plotting the great circle points.
ax.plot(longitudes, latitudes, 'r', transform=ccrs.Geodetic())
plt.plot(lon1, lat1, 'bo', markersize=7, transform=ccrs.Geodetic())
plt.plot(lon2, lat2, 'ro', markersize=7, transform=ccrs.Geodetic())

plt.show()

解决方案 #

问题根源是自定义的great_circle_points函数中，大圆坐标的计算逻辑有误，导致生成的路径点偏离实际大圆轨迹。以下提供两种修复方式：

方法一：使用Cartopy内置功能绘制大圆路径（推荐） #

Cartopy已内置处理大圆路径的功能，无需手动计算坐标，直接通过ax.plot并指定transform=ccrs.Geodetic()即可自动生成两点间的准确大圆路径，代码更简洁可靠。

修改后的代码：

import matplotlib.pyplot as plt
import cartopy.crs as ccrs

lat1, lon1 = -33.9715, 18.6021  # Cape Town
lat2, lon2 = -42.8364, 147.5075  # Hobart

plt.figure(figsize=[10, 10])
ax = plt.axes(projection=ccrs.AzimuthalEquidistant(central_latitude=90, central_longitude=0))
ax.set_global()
ax.coastlines()

# 直接绘制两点间的大圆路径
ax.plot([lon1, lon2], [lat1, lat2], 'r', transform=ccrs.Geodetic())
# 绘制端点标记
ax.plot(lon1, lat1, 'bo', markersize=7, transform=ccrs.Geodetic())
ax.plot(lon2, lat2, 'ro', markersize=7, transform=ccrs.Geodetic())

plt.show()

方法二：修复自定义的大圆坐标计算函数 #

若需自行实现大圆坐标计算，需修正函数中的球面插值公式，使用正确的球面线性插值（SLERP）逻辑：

修正后的完整代码：

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

def great_circle_points(lat1, lon1, lat2, lon2, num_points):
    # 将坐标转换为弧度
    lat1, lon1, lat2, lon2 = np.radians([lat1, lon1, lat2, lon2])
    fractions = np.linspace(0, 1, num_points)
    
    # 计算两点间的圆心角
    cos_delta = np.sin(lat1)*np.sin(lat2) + np.cos(lat1)*np.cos(lat2)*np.cos(lon2 - lon1)
    delta = np.arccos(np.clip(cos_delta, -1, 1))  # 限制范围避免数值误差
    
    # 球面线性插值计算路径点
    sin_delta = np.sin(delta)
    sin_1_minus_t_delta = np.sin((1 - fractions)*delta)
    sin_t_delta = np.sin(fractions*delta)
    
    lat = np.arctan2(
        sin_1_minus_t_delta*np.sin(lat1) + sin_t_delta*np.sin(lat2),
        sin_1_minus_t_delta*np.cos(lat1) + sin_t_delta*np.cos(lat2)*np.cos(lon2 - lon1)
    )
    lon = lon1 + np.arctan2(
        sin_t_delta*np.cos(lat2)*np.sin(lon2 - lon1),
        sin_1_minus_t_delta*np.cos(lat1) + sin_t_delta*np.cos(lat2)*np.cos(lon2 - lon1)
    )
    
    # 转换回角度坐标
    lat = np.degrees(lat)
    lon = np.degrees(lon)
    
    return lat, lon

lat1, lon1 = -33.9715, 18.6021  # Cape Town
lat2, lon2 = -42.8364, 147.5075  # Hobart
num_points = 1001

latitudes, longitudes = great_circle_points(lat1, lon1, lat2, lon2, num_points)

plt.figure(figsize=[10, 10])
ax = plt.axes(projection=ccrs.AzimuthalEquidistant(central_latitude=90, central_longitude=0))
ax.set_global()
ax.coastlines()

ax.plot(longitudes, latitudes, 'r', transform=ccrs.Geodetic())
ax.plot(lon1, lat1, 'bo', markersize=7, transform=ccrs.Geodetic())
ax.plot(lon2, lat2, 'ro', markersize=7, transform=ccrs.Geodetic())

plt.show()

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