1

我正在尝试创建一个夜灯的 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)
4

2 回答 2

1

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

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

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

persp(volcano)

在此处输入图像描述

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

在此处输入图像描述

于 2020-07-12T05:47:06.210 回答
1

你可以做

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

persp(r)

或者如果你想与它互动

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

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

在此处输入图像描述

于 2020-07-12T06:01:32.290 回答