Python极坐标散点数据绘制3D曲面plot_surface格式适配问题

核心问题 #

plot_surface方法要求输入的X、Y、Z三个参数均为同维度的二维矩阵，你当前的x、y、z都是一维散点序列，直接对一维x、y调用np.meshgrid生成的二维网格和一维z值无法对应，因此无法满足输入要求。

两种可行解决思路 #

• 方案1：利用你现有极坐标规则采样的特性直接构建二维矩阵，无需额外插值
  你给出的radius和degree是规则采样：相同半径对应固定8个角度采样点，只要先提取唯一的半径、角度列表，构建极坐标二维网格，再将对应的z值按位置填充为二维矩阵，最后转笛卡尔坐标即可直接调用plot_surface。
• 方案2：针对不规则散点通用插值方案
  如果是无规则的散点数据，可以用scipy.interpolate.griddata将一维x/y/z插值到自定义的规则二维网格上，再生成曲面。

可直接运行的修改后代码 #

from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
import numpy as np
import math
 
def polar_to_cartesian(ih, r):
    # 直接支持numpy数组运算，无需循环
    x = ih * np.cos(np.deg2rad(r))
    y = ih * np.sin(np.deg2rad(r))
    return x, y

# 原始数据转numpy数组
radius = np.array([0, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.751, 0.751, 0.751, 0.751, 0.751, 0.751, 0.751, 0.751, 1.251, 1.251, 1.251, 1.251, 1.251, 1.251, 1.251, 1.251, 1.501, 1.501, 1.501, 1.501, 1.501, 1.501, 1.501, 1.501, 2.002, 2.002, 2.002, 2.002, 2.002, 2.002, 2.002, 2.002, 2.252, 2.252, 2.252, 2.252, 2.502, 2.502, 2.502, 2.502])
degree = np.array([0, 22.5, 157.5, 202.5, 337.5, 67.5, 112.5, 247.5, 292.5, 22.5, 157.5, 202.5, 337.5, 67.5, 112.5, 247.5, 292.5, 22.5, 157.5, 202.5, 337.5, 67.5, 112.5, 247.5, 292.5, 22.5, 157.5, 202.5, 337.5, 67.5, 112.5, 247.5, 292.5, 22.5, 157.5, 202.5, 337.5, 67.5, 112.5, 247.5, 292.5, 22.5, 157.5, 202.5, 337.5, 22.5, 157.5, 202.5, 337.5])
z_hor = np.array([0, -1.23, -2.89, 2.62, -1.305, -5.75, -8.975, -1.35, -1.975, -2.675, -6.925, 2.81, -3.47, -6.49, -13.1, -0.18, -3.4, -6.735, -11.02, 2.705, -6.31, -9.96, -18.535, -0.785, -3.735, -11.64, -17.855, 0.47, -12.46, -15.195, -28.37, -2.1, -7.49, -12.46, -18.54, 2.85, -11.97, -15.08, -10.215, -10.28, -19.045, -4.26, -15.405, -6.24, -1.965, -14.78, -19.235, 3.36, -13.695])
z_ver = np.array([0, -3.075, -1.715, 3.995, 0.22, -4.41, -3.855, 3, 2.42, -6.015, -1.435, 2.64, 1.31, -8.385, -5.055, 0.825, 3.465, -7.44, -1.76, 2.695, 2.81, -12.61, -7.97, 0.575, 6.14, -15.2, -4.595, -1.23, 0.485, -20.855, -13.245, -3.89, 4.805, -12.92, -0.18, 2.345, 5.41, -18.59, -17.825, -10.485, 6.88, -1.665, -1.9, -2.14, -0.51, -14.245, 1.78, 9.575, 8.33]) 

# 提取唯一的半径和角度，按顺序排序
r_unique = np.unique(radius)
theta_unique = np.unique(degree)

# 构建二维Z矩阵，这里以z_hor为例，z_ver同理
Z = np.zeros((len(r_unique), len(theta_unique)))
for i, r in enumerate(r_unique):
    for j, theta in enumerate(theta_unique):
        # 找到对应r和theta的z值
        idx = np.where((radius == r) & (degree == theta))[0]
        if len(idx) > 0:
            Z[i,j] = z_hor[idx[0]]

# 极坐标网格转笛卡尔坐标
R, Theta = np.meshgrid(r_unique, theta_unique, indexing='ij')
X, Y = polar_to_cartesian(R, Theta)

# 绘制曲面
fig = plt.figure(figsize=(10,8))
ax1 = fig.add_subplot(111, projection='3d')

# 绘制曲面，设置透明度方便观测
surf = ax1.plot_surface(X, Y, Z, cmap='viridis', alpha=0.7)
# 可同时保留原来的散点做对照
ax1.scatter(radius*np.cos(np.deg2rad(degree)), radius*np.sin(np.deg2rad(degree)), z_hor, c='red', s=20)

ax1.set_xlabel('x [mm]')
ax1.set_ylabel('y [mm]')
ax1.set_zlabel('z[µm]')
# 添加色标
fig.colorbar(surf, shrink=0.5, aspect=5)

plt.show()

补充说明 #

如果要绘制z_ver的曲面，只需要把代码中填充Z矩阵时的z_hor替换为z_ver即可。如果后续采样为不规则散点，只需引入scipy.interpolate.griddata，将一维x/y/z插值到自定义的X/Y网格上得到二维Z矩阵即可调用plot_surface。

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