R中的表面栅格图答案

【问题标题】：Surface raster map in RR中的表面栅格图
【发布时间】：2020-07-12 04:00:43
【问题描述】：

我正在尝试创建一个夜灯的 3D 光栅图，其值表示为高程（表面图，尖峰表示值）。我已经裁剪了栅格，到目前为止，我的代码以 2D 的形式投影地图，如下所示，但不确定如何添加 3D 或应该使用哪个函数。

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()

更新！

我按照建议使用 plotly 的说明进行操作，但得到一个空矩阵。

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)

我测试了 persp(raster_new)，它给了我一个我想要的类似结果，但没有 3D 运动，小而全黑（所以没有用）。这是我的光栅的尺寸：

persp(raster_new)

我不确定用 plotly 处理矩阵时出了什么问题。在我的栅格信息下方

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)

【问题讨论】：

如果您包含一个简单的reproducible example，其中包含可用于测试和验证可能解决方案的示例输入和所需输出，则更容易为您提供帮助。

标签： r dictionary surface 3d-mapping

【解决方案1】：

Base R 提供persp() 用于矩阵的 3D 可视化，plotly 提供了很好的交互式 3D 功能。

这意味着您需要将类矩阵的对象传递给这些函数。我相信您可以将RasterLayer 从raster 包中强制转换为带有as.matrix(my_raster) 的矩阵。使用class(my_converted_raster) 检查它确实是一个矩阵。

使用内置的volcano 数据（已经是矩阵）：

persp(volcano)

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

【讨论】：

感谢@Rich PaulooIt，该绘图适用于火山数据，但对于我的矩阵栅格，它只显示一个空图。查看更新。

【解决方案2】：

你可以的

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

persp(r)

或者如果你想与之互动

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

这会在您的计算机上打开一个 rgl 设备：

【讨论】：