Surface raster map in R

Question

I am trying to create a 3D raster map of nightlights with the values represented as an elevation (surface map with spikes representing the value). I already crop the raster and so far my code to project the map, in 2D looks, like below but not sure how to add the 3D or which function should I use.

ras_df <- as.data.frame(crop_raster, xy=TRUE, na.rm=TRUE)
colnames(ras_df) <- c("x", "y", "value")
ras_df$log_value_plus1 <- log(ras_df$value+1)

fig <- ggplot() +
  geom_raster(data=ras_df, aes(x=x,y=y, fill=log_value_plus1), alpha=1) +
 
  geom_polygon(data=shapfile_fortified, aes(x=long, y=lat, group=group), 
               fill=NA, color="grey60", size=0.25) +
   coord_quickmap()
  

UPDATE!

I follow the instructions suggested using plotly but I get an empty matrix.

raster_new<-raster::as.matrix(crop_raster)
class(raster_new)
#matrix
plot_ly(z = ~raster_new) %>% add_surface()

#same empty matrix using
plot_ly(z = raster_new,  type="surface",showscale=FALSE)

I tested persp(raster_new) and it gives me a similar result I would like to have but with no 3D movement, small and all dark (so no useful). This are the dimensions of my raster:

persp(raster_new)

I am no sure what is going wrong when the matrix is processed with plotly. Below my raster info

class      : RasterLayer 
dimensions : 4352, 4959, 21581568  (nrow, ncol, ncell)
resolution : 0.004166667, 0.004166667  (x, y)
extent     : -8.672916, 11.98958, 18.96042, 37.09375  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : light_raster_2016_.tif 
names      : light_raster_ALG_2016_ 
values     : 0, 21980.4  (min, max)

Answer 1

Base R provides persp() for 3D visualization of matrices, and plotly offers a nice interactive 3D capabilities.

This means you'll need to pass an object of class matrix to these functions. I believe you can coerce a RasterLayer from the raster package into a matrix with as.matrix(my_raster). Use class(my_converted_raster) to check that it indeed is a matrix.

Using the built-in volcano data (already a matrix):

persp(volcano)

library(plotly)
plot_ly(z = ~volcano) %>% add_surface()

Answer 2

You can do

library(raster)
r <- raster(volcano)
extent(r) <- c(0, 610, 0, 870)

persp(r)

Or if you want to interact with it

library(rasterVis)
library(rgl)
plot3D(r, col = rainbow)

Which opens an rgl device on your computer: