rrgl

How can I plot contour lines on a sphere?

Suppose I have an rgl spherical object and I want to plot contours on the surface of this sphere. I figure that the contourLines3d function from the rgl R package is my best bet, so I'm going to adapt their MWE for us here:

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 # radius of Earth in 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(ids, func)

where I need to specify func properly to plot the contours. This is where I'm struggling. Suppose I have a pdf function (dSphereFunc) that computes the density of any point on the sphere. How can I plot contours using my dSphereFunc? I'm sure there's a relatively simple way to do this, but I'm really having troubles here. Thank you for your help.

Update: As requested in the comments, I have supplied the density function that I'm trying to evaluate:

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) # formatting from matrix to constant
    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)
}

In the above function, we are treating mu (a dx1 vector) and V (a dxd matrix) as fixed (i.e., already estimated), where d=3; additionally, y is a dx1 vector which is the spherical coordinates x, y, and z.

To have a MWE, here are some values for mu and Vinv:

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))
test=c(1,3,4)
test=test/norm(test,type="2")
dSphereFunc(test,mu,Vinv) # 1.156622

I think my trouble is that dSphereFunc takes three arguments, the latter two being fixed. So, I tried using the following function which only requires one argument, xyz which is a dx1 vector for the spherical coordinates x, y, and z:

densFunc = function(xyz){apply(xyz, 1, function(point) dSphereFunc(point, mu, Vinv))}

However, this is still not working.

Solution

  • The code you were using was fine, the only problem is with the plotting: the contours all appear on the back of the sphere, so you can't see them. Setting alpha = 0.5 when you draw the sphere makes them available, or you can rotate the sphere until they show up. Here the lines are showing through from the back of the sphere:

    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 # radius of Earth in km
x <- r*cos(lat)*cos(long)
y <- r*cos(lat)*sin(long)
z <- r*sin(lat)

library(rgl)
open3d()
#> glX 
#>   1
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, alpha = 0.5)

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) # formatting from matrix to constant
  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))
test=c(1,3,4)
test=test/norm(test,type="2")
dSphereFunc(test,mu,Vinv) 
#>          [,1]
#> [1,] 1.156622

densFunc = function(xyz){apply(xyz, 1, function(point) dSphereFunc(point, mu, Vinv))}


contourLines3d(ids, densFunc)

    Created on 2024-05-09 with reprex v2.1.0