How to project a 2D image as a surface in a 3D scatter plot (in R)?

I'm working with 3 dimensional coordinates data, which i'm plotting in a scatterplot, i have ~30.000 datapoints, and i've included the first 10 here so that you can reproduce it

library(rgl)
library(plot3D)
library(car)

df <- data.frame(meanX = c(147.34694,
                       173.89244,
                       135.73004,
                       121.93766,
                       109.72152,
                       92.53709,
                       165.46588,
                       169.77744,
                       127.01796,
                       99.34347),
             meanY = c(140.40816,
                       110.99128,
                       134.56023,
                       164.18703,
                       166.04051,
                       155.97329,
                       105.29377,
                       104.42683,
                       130.17066,
                       155.99696),
             avgDist = c(40.788118,
                         12.957329,
                         14.24348,
                         39.10424,
                         34.694258,
                         25.532335,
                         21.491695,
                         23.528944,
                         9.309201,
                         31.916879))

I've been using the scatter3d function to plot this

scatter3d(x = df$meanX, y = df$meanY, z = df$avgDist, surface = FALSE)

Now my "problem", is that I would like to have a 2d surface with an external image file overlayed onto it at z=0, and as a bonus, if i could project a heatmap/contours from the scatterplot data (meanX and meanY used for the contours) over that image as well, that would be great.

This is the image i'd like to have draped at z = 0:

https://i.sstatic.net/BOnnv.png

That image was made with this ggplot:

map.colors <- colorRampPalette(c("green","yellow","red"))

densityPlot <- ggplot(direData, aes(x = meanX, y = ,meanY)) + 
  stat_density2d(geom="tile", aes(fill=..density.., alpha=sqrt(sqrt(..density..))), contour=FALSE, n=100) +
  scale_alpha(range = c(0, 1.0)) + scale_fill_gradientn(colours = map.colors(5)) + 
  xlim(70,185) + ylim(70,185)

minimap <- readPNG('~/yasp/minimap.png')

densityPlot + annotation_raster(minimap, ymin = 70 ,ymax=185 ,xmin = 70,xmax = 185) + 
  stat_density2d(geom="tile", aes(fill=..density.., alpha=10*sqrt(..density..)), contour=FALSE, n=100)

Is there any way to do this? I've googled quite a bit for a solution but found no real way of doing this. I don't mind creating the image first in ggplot2 with the heatmap, saving that, and then using that as input for the surface draping, but it would of course be quite cool if it could all be done in one call to plot.

Solution

(2nd Edit) I try to write something better code and confirm two xy-coordinates are the same. ggplot2 theme with no axes or grid help me to plot only the panel region.

library(rgl); library(grid); library(gtable)

df <- data.frame(meanX = c(147.34694, 173.89244, 135.73004, 121.93766,
                           109.72152,  92.53709, 165.46588, 169.77744,
                           127.01796,  99.34347),
                 meanY = c(140.40816, 110.99128, 134.56023, 164.18703,
                           166.04051, 155.97329, 105.29377, 104.42683,
                           130.17066, 155.99696),
                 avgDist = c(40.788118, 12.957329, 14.24348, 39.10424,
                             34.694258, 25.532335, 21.491695,23.528944,
                             9.309201,  31.916879))

map.colors <- colorRampPalette(c("green","yellow","red"))

# set scale_*_continuous() to plot only the panel region. limits mean xlim (or ylim)
# change "tile" into "raster" because of making noise lines on my screen
densityPlot <- ggplot(df[,1:2], aes(x = meanX, y = ,meanY)) + 
  stat_density2d(geom="raster", aes(fill=..density.., alpha=sqrt(sqrt(..density..))), contour=FALSE, n=100) +
  scale_alpha(range = c(0, 1.0)) + scale_fill_gradientn(colours = map.colors(5)) + 
  scale_x_continuous(limits=c(70,185), expand = c(0,0)) + scale_y_continuous(limits=c(70,185), expand = c(0,0)) +
  geom_point(size=4)               # to test XY-coordinate (black points on the map)

open3d()
plot3d( df, type="s", radius=1, col="red", axes=F, 
        xlim = c(70,185), ylim = c(70,185),
        expand = 1 )
plot3d( df, type="h", col="blue", add=T )  # to test XY-coordinate (line segments from z = 0)
axes3d(c("x","y","z") )
show2d({                  # show2d uses 2D plot function's output as a texture on a box.
  grid.draw(gtable_filter(ggplotGrob(densityPlot), "panel"))
},
expand = 1 , texmipmap = F )   # texmipmap = F makes tone clear (not essential)

# I think this is clearly better than using a intermediate file,
# so I deleted related code. Thanks Mike !

plot