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.

Please hover over the R tag - it asks for a a) minimal & b) reproducible example, with which one can recreate your current state & problem using copy,paste,run. Maybe you want something like library(raster);library(plot3D);library(rgl);slogo <- stack(system.file("external/rlogo.grd", package="raster"));f <- function(x, offset=0, add=FALSE) { x<-as.matrix(x);persp3D(z=x+offset, zlim=c(0,3000), add=add, phi=50, theta=30);persp3d(z=x+offset, zlim=c(0,3000), col="lightblue", add=add, alpha = .75)};f(slogo[[1]],0*1000,F);f(slogo[[3]],1*1000,T);f(slogo[[3]],2*1000,T)? — lukeA
@lukeA I still have some problems. Why the 3D plot is changing the original shape of the polygon? (I have attached the code lines for a reproducible example above). Just plotted the "r1" raster with a simple "plot" function and yours, and it seems that something is happening with the polygon. — ADM

lukeA lukeA · Accepted Answer · 2017-04-16T12:13:13

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)

3D raster stack plot

1 Answers