如何用R语言rgl包在球体表面绘制自定义密度轮廓线？

核心问题分析 #

你的dSphereFunc要求输入是单位球面上的点，但contourLines3d传入的是你创建的地球球面实际坐标（半径6378.1），未做归一化；同时需要确保func参数严格符合contourLines3d的输入输出要求：接收N×3的坐标矩阵，返回长度为N的数值向量。

修正后的代码实现 #

# 自定义密度函数（保持原逻辑不变）
dSphereFunc = function(y, mu, Vinv){
    d = length(mu)
    ytVinvy = t(y)%*%Vinv%*%y
    ytmu = t(y)%*%mu
    result = ((2*pi)^(-(d-1)/2))/(ytVinvy^(d/2))*exp((ytmu^2/ytVinvy - t(mu)%*%mu)/2)*mdmo(d, ytmu/(ytVinvy^(1/2)))
    return(result)
}

mdmo = function(d, t){
    t = c(t) # 格式化为常量
    integrand = function(x) x^(d-1)*exp(-(x-t)^2/2)
    integral = integrate(integrand, lower = 0, upper = Inf)
    return((2*pi)^(-1/2)*integral$value)
}

# 固定参数
mu=c(2,2,2)
Vinv=solve(matrix(c(1.57,-0.08,-0.5,-0.08,0.74,0.34,-0.5,0.34,1.16),nrow=3))

# 修正后的密度函数：处理坐标归一化，符合contourLines3d要求
densFunc = function(xyz){
    # 对每个点归一化到单位球
    unit_points = t(apply(xyz, 1, function(point) point / norm(point, type="2")))
    # 计算每个单位点的密度
    apply(unit_points, 1, function(point) dSphereFunc(point, mu, Vinv))
}

# 生成球面坐标（原代码不变）
lat <- matrix(seq(90, -90, len = 50)*pi/180, 50, 50, byrow = TRUE)
long <- matrix(seq(-180, 180, len = 50)*pi/180, 50, 50)
r <- 6378.1 # 地球半径（km）
x <- r*cos(lat)*cos(long)
y <- r*cos(lat)*sin(long)
z <- r*sin(lat)

# 绘制球面并添加轮廓线
open3d()
ids <- persp3d(x, y, z, col = "white",
        specular = "black", axes = FALSE, box = FALSE, xlab = "", ylab = "", zlab = "",
        normal_x = x, normal_y = y, normal_z = z, polygon_offset = 1)

# 调用contourLines3d，使用修正后的densFunc
contourLines3d(ids, densFunc, levels = seq(0.1, 1.2, by = 0.2)) # 可根据密度范围调整levels

关键修正点 #

• 坐标归一化：在densFunc中先将传入的球面坐标转换为单位向量，匹配dSphereFunc的输入要求。
• 函数格式兼容：确保返回的是长度与输入行数一致的数值向量，apply默认按行处理后返回向量，符合contourLines3d的要求。
• 自定义轮廓层级：通过levels参数指定要绘制的轮廓线对应的密度值，可根据实际计算的密度范围调整。

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