Python极坐标插值后图面中心存在缺失点解决方法

缺失的核心原因 #

• 核心错误根源是将极坐标下不同物理意义的两个维度（角度、径向距离）等价为平面直角坐标输入griddata做三角剖分。两个维度量纲不统一、几何意义错配，会导致Delaunay算法计算出的凸包形状，和实际测量的圆形区域完全不匹配。
• 从测点分布来看，除圆心孤立测点外，其余所有测点全部分布在r=75、r=144两个同心圆上。在（角度, 径向距离）组成的二维平面上，凸包的下边界是连接t=π、r=75位置与圆心（t=0, r=0）的两条直线，这会导致t接近±π的内侧低径向距离区域完全落在凸包范围外。而griddata的线性、三次插值方法，对凸包外的点会直接返回空值（nan），最终在图上形成连片的中心空洞。
• 原代码中仅手动添加角度周期重复点的逻辑，既无法补上低径向距离区域的凸包覆盖缺口，也无法缓解r=0位置的坐标奇异带来的剖分误差。

具体修正方案 #

• 切换插值平面：放弃在极坐标平面开展插值的思路，直接采用原始测点的笛卡尔坐标（x, y）作为插值输入。这种坐标下的Delaunay剖分，凸包就是测点所在的半径144的圆形测量区域，完全覆盖目标范围，从根源上解决几何错配问题。
• 网格适配：生成目标极坐标网格后，将每个网格点转换为对应的笛卡尔坐标值，作为插值的目标位置，不需要在极坐标空间内执行剖分计算。
• 空值兜底：线性插值仅会对凸包边缘的极少数点返回空值，使用最近邻插值的结果对所有空值位置做填补，保证整个圆形区域无缺失像素。
• 绘图对齐：Matplotlib的极坐标角度默认范围为0到2π，绘图前将负角度值转换为对应的正角度值，避免填充区域、测点标注的位置偏移。

修正后的可运行代码 #

import numpy as np
from scipy.interpolate import griddata
import matplotlib.pyplot as plt

# 原始测量测点数据
x = np.array([-144.00, -101.80, -101.80, -75.00, -53.00, -53.00, 0.00, 0.00, 0.00, 0.00, 0.00, 53.00, 53.00, 75.00, 101.80, 101.80, 144.00])
y = np.array([0.00, 101.80, -101.80, 0.00, 53.00, -53.00, 144.00, 75.00, 0.00, -75.00, -144.00, 53.00, -53.00, 0.00, 101.80, -101.80, 0.00])
z = np.array([148.3861807, 148.9051447, 148.0415147, 147.9976293, 147.98485, 147.9579673, 148.89261, 148.0217707, 147.9312247, 147.7952, 147.3225247, 148.2489567, 148.2120013, 148.3169953, 149.578092, 147.9356893, 148.556672])

# 生成极坐标下的目标网格
r_grid, t_grid = np.meshgrid(np.linspace(0, 144, 50), np.linspace(-np.pi, np.pi, 50))
# 将极坐标网格转换为笛卡尔坐标，作为插值目标位置
x_grid = r_grid * np.cos(t_grid)
y_grid = r_grid * np.sin(t_grid)

# 笛卡尔平面上执行插值
# 线性插值获得主体平滑结果
z_linear = griddata((x, y), z, (x_grid, y_grid), method='linear')
# 最近邻插值结果兜底填补空值
z_nearest = griddata((x, y), z, (x_grid, y_grid), method='nearest')
z_final = np.where(np.isnan(z_linear), z_nearest, z_linear)

# 绘制极坐标填充图
fig, ax = plt.subplots(figsize=(10,10), subplot_kw=dict(projection='polar'))
# 角度转换为0-2pi范围，匹配Matplotlib极坐标规则
t_plot = np.mod(t_grid, 2 * np.pi)
ax.contourf(t_plot, r_grid, z_final, 50, cmap='jet')

# 标注原始测点的数值
for xi, yi, zi in zip(x, y, z):
    ti = np.arctan2(yi, xi)
    ri = np.sqrt(xi**2 + yi**2)
    ti_plot = np.mod(ti, 2 * np.pi)
    plt.text(ti_plot, ri, "{:.2f}".format(zi), ha="center", va="center", color="k", fontsize=8)

plt.show()

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