如何基于非完整网格绘制Surface Plot？

你遇到的核心问题是：plot3D包的mesh()、slice3D()、isosurf3D()这些函数都要求输入规则的网格数据——也就是x和y的每一个组合都必须有对应的z值，但你的数据是缺失大量组合的非规则（稀疏）数据，所以这些函数无法正常运行。

下面是具体的解决方案，我们可以通过插值补全网格来生成规则数据，再绘制Surface Plot：

步骤1：整理数据并安装必要的包 #

首先把你的数据整理成数据框，然后安装并加载用于插值的akima包（也可以用fields等其他插值包）：

# 整理数据
df <- data.frame(
  x = c(10L, 20L, 30L, 40L, 50L, 60L, 70L, 80L, 90L, 100L, 30L, 40L, 50L, 60L, 70L, 80L, 90L, 100L, 50L, 60L, 70L, 80L, 90L, 100L, 70L, 80L, 90L, 100L, 90L, 100L),
  y = c(10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 20L, 20L, 20L, 20L, 20L, 20L, 20L, 20L, 30L, 30L, 30L, 30L, 30L, 30L, 40L, 40L, 40L, 40L, 50L, 50L),
  z = c(6.093955007, 44.329214443, 149.103755156, 351.517349974, 726.51174655, 1191.039562104, 1980.245204702, 2783.308022984, 6974.563519067, 5149.396230019, 142.236259009, 321.170609648, 684.959503897, 1121.475597135, 1878.334840961, 2683.116309688, 4159.60732066, 5294.774284119, 687.430547359, 1119.765405426, 1876.57337196, 2685.951176024, 3945.696884503, 5152.986796572, 1870.78724464, 2677.744176903, 3951.928931107, 5160.295960254, 3957.503273558, 5147.237754092)
)

# 安装并加载插值包
install.packages("akima")
library(akima)
library(plot3D)

步骤2：插值生成规则网格 #

用akima::interp()函数对稀疏数据进行插值，生成完整的x-y网格和对应的z值：

# 生成规则网格的x和y序列（可以调整n值控制精度）
x_seq <- seq(min(df$x), max(df$x), length.out = 20)
y_seq <- seq(min(df$y), max(df$y), length.out = 20)

# 插值
interp_result <- interp(x = df$x, y = df$y, z = df$z, 
                        xo = x_seq, yo = y_seq, 
                        linear = TRUE) # 线性插值，也可以用spline插值

步骤3：绘制Surface Plot #

现在有了规则网格数据，就可以用plot3D::persp3D()绘制Surface Plot了：

# 绘制Surface Plot
persp3D(x = interp_result$x, y = interp_result$y, z = interp_result$z,
        xlab = "x", ylab = "y", zlab = "z",
        col = "lightblue", border = "gray",
        theta = 35, phi = 20, # 调整视角
        main = "Surface Plot with Interpolated Data")

补充：修复slice3D和isosurf3D的运行问题 #

如果要运行你原来的slice3D()和isosurf3D()，同样需要用插值后的规则网格数据：

# 基于插值后的网格计算p值
M <- mesh(interp_result$x, interp_result$y)
R <- with(M, sqrt(x^2 + y^2 + interp_result$z[cbind(M$ind[1,], M$ind[2,])]))
p <- sin(2*R)/(R+1e-3)

# 转换为矩阵格式
p_matrix <- matrix(p, nrow = length(x_seq), ncol = length(y_seq))

# 运行slice3D
slice3D(x = interp_result$x, y = interp_result$y, z = interp_result$z,
        colvar = p_matrix, xs = 50, ys = c(10,30,50), zs = NULL)

# 运行isosurf3D
isosurf3D(x = interp_result$x, y = interp_result$y, z = interp_result$z,
          colvar = p_matrix, level = 0, col = "red")

注意事项 #

• 插值方法的选择：linear线性插值适合数据变化平缓的情况，如果你的数据有明显的非线性趋势，可以尝试spline插值（把interp()里的linear设为FALSE）。
• 网格精度：length.out的值越大，网格越精细，但计算量也会增加，你可以根据需求调整。

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