3D raster stack plot

Question

I have a raster stack with 2 layers which correspond to 2 elevations. Every layer is a raster with Lon and Lat coordinates with value 1 (all missing values are NA), forming polygons. RasterStack example

r1<-raster(xmn = 2.5, xmx = 3.3, ymn = 42, ymx = 42.5, nrows = 28, ncols = 33)
r2<-raster(xmn = 2.5, xmx = 3.3, ymn = 42, ymx = 42.5, nrows = 28, ncols = 33)

a<-c(421, 422, 424, 453, 454, 455, 456, 457, 485, 486, 487, 488, 489, 513, 514, 515, 516, 517, 518, 519, 546, 547, 548, 549, 550, 579, 580, 581, 582, 583, 613, 614, 615, 646, 647, 648, 649, 680, 681, 682)
r1[a]<-1
b<-c(514, 515, 516, 547, 548, 549, 550, 613, 614, 615, 647, 648, 649)
r2[b]<-1

st<-stack(r1,r2)

Is there any way to make a 3D plot showing every polygon in a different elevation depending on the layer? (x=LON, Y=LAT, Z= nlayer)

I have achieve a scatter3D plot for layers 1 and 2 Scatter3D example, but I would like to show them as polygons, to see how areas overlap within the entire vertical.

Answer 1

You could try

library(raster)
library(plot3D)
f <- function(x, offset=0, add=FALSE) { 
  poly <- rasterToPolygons(x, dissolve = T)
  coords <- poly@polygons[[1]]@Polygons[[1]]@coords
  polygon3D(
    x=coords[,1], 
    y=coords[,2], 
    z=rep(offset,nrow(coords)), 
    xlim=lims[1:2],
    ylim=lims[3:4],
    zlim=lims[5:6], 
    add=add, 
    col="red",
    theta = 0
  )
}
lims <- c(as(extent(st),"vector"), 0, 3000)
f(st[[1]],0*1000,F)
for (x in 2:nlayers(st)) 
  f(st[[x]],x*1000,T)