[R-sig-Geo] spatial data 3D bivariate distributions

[Paolo Piras]

Hi folks,
I'm trying to modify an R code for plotting x,y spatial data in 3D  bivariate distribution using R;
I modified a code from here:
http://stackoverflow.com/questions/19949435/3d-plot-of-bivariate-distribution-using-r-or-matlab

 I would like to obtain the same results coming from the code suggested by r2evans user but with distributions separated by groups, possibly with semi-transparent color. I have problems in:

1) plotting the histograms on y axis; I just cannot "translate" the procedure for x in the y axis
2) putting all histograms perfectly superimposed and not shifted one from each other

here below the modified r2evans's code with a stand alone reprodicible example.
thanks in advance for your suggestions
best
paolo



library(rgl)
bvn<-matrix(rnorm(900,0,0.1),ncol=3)
group<-factor(rep(1:3,each=100))
xmaxlist<-NULL
ymaxlist<-NULL
zmaxlist<-NULL
hxlist<-NULL
hxslist<-NULL
hylist<-NULL
hyslist<-NULL

for(i in 1: nlevels(group)){
  hxi<- hist(bvn[,2][as.numeric(group)==i], plot=FALSE)
  hxsi <- hxi$density / sum(hxi$density)
  hyi <- hist(bvn[,1][as.numeric(group)==i], plot=FALSE)
  hysi <- hyi$density / sum(hyi$density)
  
  ## [xy]max: so that there's no overlap in the adjoining corner
  xmaxi <- tail(hxi$breaks, n=1) + diff(tail(hxi$breaks, n=2))
  ymaxi <- tail(hyi$breaks, n=1) + diff(tail(hyi$breaks, n=2))
  zmaxi <- max(hxsi, hysi)
  
  hxlist<-c(hxlist,list(hxi))
  hxslist<-c(hxslist,list(hxsi))
  hylist<-c(hylist,list(hyi))
  hyslist<-c(hyslist,list(hysi))
  xmaxlist<-c(xmaxlist,list(xmaxi))
  ymaxlist<-c(ymaxlist,list(ymaxi))
  zmaxlist<-c(zmaxlist,list(zmaxi))
}
for(i in 1:length(zmaxlist)){
  
  for (ii in seq_along(hxlist[[i]]$counts)) {
    quads3d(hxlist[[i]]$breaks[ii]*c(.9,.9,.1,.1) + hxlist[[i]]$breaks[ii+1]*c(.1,.1,.9,.9),
            rep(ymaxlist[[i]], 4),
            hxslist[[i]][ii]*c(0,1,1,0), color=i,alpha=0.5,add=T)
  }
  
}

plot3d(bvn[,2], bvn[,1], 0, zlim=c(0, max(unlist(zmaxlist))), radius=0.02,axes=FALSE,col=as.numeric(group),type="s",add=T)