Filling a region on a sphere with rgl

Here is the spherical Voronoï tessellation of the world capitals:

I have the coordinates of the points defining the boundary of each country. I'd like to fill the countries with a color. The motivation is to plot the Voronoï edges with the same color as the seas/oceans so that we will not see them outside the countries. In other words, I'm looking for a function similar to the polygon function but for spherical polygons.

Edit

Well, I get it now. But for some reason, the cyan color of the water and the cyan color of the Voronoï edges are different:

Here is how to fill the countries.

library(maps)
data(worldMapEnv)

# make a PNG image of the world in 2D with filled countries
world <- map("world", plot = FALSE, fill = TRUE, wrap = c(0, 360))

png("world2D.png", width = 1024L, height = 1024L)
opar <- par(mar = c(0, 0, 0, 0))
plot(world$x, world$y, type = "n", xlab = NA, ylab = NA, axes = FALSE,
     xaxs = "i", yaxs = "i", xlim = c(0, 360), ylim = c(-90, 90))
polygon(world$x, world$y, col="orangered", lwd = 2)
par(opar)
dev.off()

# use spherical coordinates to make a mesh of the unit sphere
library(cgalMeshes)
sphericalCoordinates <- function(θ, ϕ){
  x <- cos(θ) * sin(ϕ)
  y <- sin(θ) * sin(ϕ)
  z <- cos(ϕ)
  rbind(x, y, z)
}
rmesh <- parametricMesh(
  sphericalCoordinates, urange = c(0, 2*pi), vrange = c(0, pi),
  periodic = c(TRUE, TRUE), nu = 1024L, nv = 1024L
)
rmesh$normals <- rmesh$vb[-4L, ]

# get the angles θ and ϕ of the vertices of the mesh
UV <- cooltools::car2sph(t(rmesh$vb[-4L, ]))
UV <- cbind(UV[, 3L], UV[, 2L])


# now load the PNG image
library(imager)
img <- load.image("world2D.png")
# take the r, g, b channels
r <- squeeze(R(img))
g <- squeeze(G(img))
b <- squeeze(B(img))

# make interpolation functions to get the colors of the UV points
library(cooltools) # to get the `approxfun2` function
x_ <- seq(0, 2*pi, length.out = 1024L)
y_ <- seq(0, pi, length.out = 1024L)
f_r <- approxfun2(x_, y_, r)
f_g <- approxfun2(x_, y_, g)
f_b <- approxfun2(x_, y_, b)

# now, interpolate the r, g, b values
UV_r <- f_r(UV[, 1L], UV[, 2L])
UV_g <- f_g(UV[, 1L], UV[, 2L])
UV_b <- f_b(UV[, 1L], UV[, 2L])

# convert rgb to hex codes
clrs <- rgb(UV_r, UV_g, UV_b)
clrs[clrs == "#FFFFFF"] <- "cyan" # replace white with cyan


# assign the colors to the vertices of the mesh
# I don't know why, but one has to reverse the colors
rmesh$material <- list(color = rev(clrs))


# plot
library(rgl)
open3d(windowRect = 50 + c(0, 0, 512, 512), zoom = 0.7)
shade3d(rmesh, meshColor = "vertices")
snapshot3d("world3D.png", webshot = FALSE)