如何在Cartopy中为经纬度散点拟合圆（理想为椭圆）？

问题说明 #

手上有一组高纬度地区的经纬度散点数据，需要拟合最优椭圆或圆。直接使用常规平面拟合工具包得到的结果明显不合理（椭圆参数出现NaN、圆半径异常），核心原因是经纬度属于球面坐标系，不能直接套用平面拟合逻辑，必须结合地球球面特性处理。

经纬度坐标数据 #

import numpy as np
coords = np.array([
    [-153.1906979 ,   62.01707771],
    [  13.05660412,   63.15537447],
    [-175.82610203,   67.11698477],
    [ -10.31730643,   61.74562855],
    [ 168.02402748,   79.60818152],
    [ -34.46162907,   65.10894426],
    [ -57.20962503,   59.49626998],
    [ 113.70202771,   68.22239091],
    [ -80.43411993,   55.6654176 ],
    [  93.77252509,   76.19392633],
    [-104.10892084,   56.68264351],
    [  66.36158188,   67.59664968],
    [-127.75176924,   57.31577071],
    [-151.83057714,   61.64142205],
    [  17.44848859,   56.02194986],
    [-176.30087703,   66.5955554 ],
    [  -5.48747931,   61.95844561],
    [ 160.22917767,   66.07650153],
    [ -27.93440014,   67.82152994],
    [ 137.09393573,   63.71148003],
    [ -53.3290508 ,   55.79699915],
    [ 109.42329666,   75.43090294],
    [ -76.59105583,   59.18143738],
    [  89.94733587,   63.50658353],
    [-100.54585734,   55.16704225],
    [  66.15810397,   64.64851675],
    [-123.65415058,   60.14507524],
    [  41.00262656,   70.67714209],
    [-145.66917977,   68.55315102],
    [  18.34306395,   67.62222778]
])

当前Cartopy绘图代码 #

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

fig = plt.figure(figsize=(20,20))
ax = fig.add_subplot(121, projection=ccrs.NearsidePerspective(
    central_longitude=0, central_latitude=90, satellite_height=30785831
))

ax.add_feature(cfeature.NaturalEarthFeature('physical', 'ocean', '50m', facecolor='#daf7f7', alpha=0.7, zorder=0))
ax.add_feature(cfeature.NaturalEarthFeature('physical', 'land', '50m', facecolor='#ebc7a4', edgecolor='black', alpha=0.7,zorder=0))

ax.set_global()
grid = ax.gridlines(draw_labels=True)
grid.xlabel_style = {'size': 20, 'color': 'black'}
grid.ylabel_style = {'size': 20, 'color': 'black'}

ax.scatter(coords[:,0], coords[:,1], c='red', s=40, zorder=1, transform=ccrs.PlateCarree())
plt.show()

已尝试方案的问题 #

• 使用常规平面椭圆拟合方法时，椭圆短轴参数出现NaN，原因是经纬度的球面特性导致平面拟合逻辑失效
• 使用skg.nsphere_fit()得到的半径为433，未考虑单位转换（实际是球面距离的无量纲值，需转为米或角度）

可行解决方案 #

核心思路 #

经纬度是球面角度坐标，必须先转换为三维笛卡尔坐标（基于地球球体模型），再在三维空间中拟合平面/椭圆，最后将结果转回经纬度坐标用于绘图。

1. 经纬度转三维笛卡尔坐标 #

def latlon_to_cartesian(lon, lat, radius=6371000):  # 地球平均半径，单位米
    lon_rad = np.radians(lon)
    lat_rad = np.radians(lat)
    x = radius * np.cos(lat_rad) * np.cos(lon_rad)
    y = radius * np.cos(lat_rad) * np.sin(lon_rad)
    z = radius * np.sin(lat_rad)
    return np.array([x, y, z]).T

# 转换所有点
cart_coords = latlon_to_cartesian(coords[:,0], coords[:,1])

2. 拟合球面圆 #

球面圆对应三维空间中平面与地球球面的交线，步骤如下：

from sklearn.decomposition import PCA

# 用PCA拟合包含所有点的最优平面
pca = PCA(n_components=3)
pca.fit(cart_coords)
normal = pca.components_[2]  # 平面法向量（方差最小的方向）
center_cart = pca.mean_      # 平面的中心（数据的三维重心）

# 计算球面圆参数
R = 6371000  # 地球半径
d = np.abs(np.dot(normal, center_cart))  # 地心到平面的距离
spherical_radius_m = R * np.arcsin(np.sqrt((R**2 - d**2))/R)  # 球面圆的弧长半径（米）
spherical_radius_deg = np.degrees(spherical_radius_m / R)  # 转换为角度半径

3. 拟合平面椭圆（投影到球面） #

若需要拟合椭圆，先将三维点投影到拟合的平面，转为二维坐标后再拟合：

# 构建平面的正交基
u = pca.components_[0]
v = pca.components_[1]

# 将三维点投影到平面，得到二维坐标
plane_coords = np.array([
    np.dot(coord - center_cart, u) for coord in cart_coords
]), np.array([
    np.dot(coord - center_cart, v) for coord in cart_coords
])
plane_coords = np.array(plane_coords).T

# 拟合椭圆方程：Ax² + Bxy + Cy² + Dx + Ey + F = 0
def fit_ellipse(x, y):
    x = x[:, np.newaxis]
    y = y[:, np.newaxis]
    D = np.hstack((x*x, x*y, y*y, x, y, np.ones_like(x)))
    S = np.dot(D.T, D)
    C = np.zeros([6,6])
    C[0,2] = C[2,0] = 2
    C[1,1] = -1
    eigvals, eigvecs = np.linalg.eig(np.dot(np.linalg.inv(S), C))
    idx = np.nonzero(np.abs(eigvals) < 1e-8)[0]
    return eigvecs[:, idx].ravel()

ellipse_params = fit_ellipse(plane_coords[:,0], plane_coords[:,1])

4. 将拟合结果转回经纬度并绘图 #

以球面圆为例，生成圆上的点并绘图：

def cartesian_to_latlon(x, y, z, radius=6371000):
    lat = np.degrees(np.arcsin(z / radius))
    lon = np.degrees(np.arctan2(y, x))
    return lon, lat

# 生成球面圆上的点
theta = np.linspace(0, 2*np.pi, 100)
circle_plane = spherical_radius_m * np.array([np.cos(theta), np.sin(theta)]).T
circle_cart = center_cart + np.dot(circle_plane, np.array([u, v]))
circle_lon, circle_lat = cartesian_to_latlon(circle_cart[:,0], circle_cart[:,1], circle_cart[:,2])

# 绘图
fig = plt.figure(figsize=(20,20))
ax = fig.add_subplot(111, projection=ccrs.NearsidePerspective(
    central_longitude=0, central_latitude=90, satellite_height=30785831
))

ax.add_feature(cfeature.NaturalEarthFeature('physical', 'ocean', '50m', facecolor='#daf7f7', alpha=0.7, zorder=0))
ax.add_feature(cfeature.NaturalEarthFeature('physical', 'land', '50m', facecolor='#ebc7a4', edgecolor='black', alpha=0.7,zorder=0))

ax.set_global()
grid = ax.gridlines(draw_labels=True)
grid.xlabel_style = {'size': 20, 'color': 'black'}
grid.ylabel_style = {'size': 20, 'color': 'black'}

ax.scatter(coords[:,0], coords[:,1], c='red', s=40, zorder=1, transform=ccrs.PlateCarree())
ax.plot(circle_lon, circle_lat, c='blue', linewidth=3, zorder=2, transform=ccrs.PlateCarree())

plt.show()

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