在R语言中计算多边形（shapefile）的Z值

Question

我的目标是在 R 中创建一个 3D 可视化。我有一份柏林市内城区（Ortsteile）的 shapefile 文件，并想要将人口密度（每平方公里的居民数）作为 z 值进行突出显示。我已经将 shapefile 文件导入 R 中，并按以下方式对密度值（"Einwohnerd"）进行了着色：

library(rgdal)
library(sp)

berlin=readOGR(dsn="C...etc.", layer="Ortsteile")

berlin@data

col <- rainbow(length(levels(berlin@data$Name)))
spplot(berlin, "Einwohnerd", col.regions=col, main="Ortsteil Berlins", sub="Datensatz der Stadt Berlin", lwd=.8, col="black")

如何将特定多边形（城市区域）与Z值（人口/平方公里）相关联，以及如何突出显示此Z值？

希望有人能给出答案！ 最好的问候 SB

感谢您的答案，但我仍在寻找最佳方法使用密度作为Z值，以便我可以创建3D模型。 我发现无法使用形状的多边形，但可以光栅化多边形并使用矩阵进行不同的透视和旋转。

这是代码，但最终的3D可视化效果不够清晰和好。也许更好的方法是以另一种方式计算Z值，使第一组值不会如此高，或者使用多边形的中心，然后在z方向上绘制柱：

library(rgdal)
library(sp)

setwd("C:\\...")
berlin=readOGR(dsn="C:\\...\\Ortsteile", layer="Ortsteile") 

col <- rainbow(length(levels(berlin@data$Name)))  
spplot(berlin, "Einwohnerd", col.regions=col, main="Ortsteil Berlins",                 
sub="Datensatz    der Stadt Berlin", lwd=.8, col="black")

library(raster)

raster <- raster(nrows=100, ncols=200, extent(berlin)) 

test <- rasterize(berlin, raster, field="Einwohnerd")

persp(test, theta = 40, phi = 40, col = "gold", border = NA, shade = 0.5)  

for(i in seq(0,90,10)){     
persp(test, theta = 40, phi = i, col = "gold", border = NA, shade = 0.5)
}

library(rgl)         
library(colorRamps)
mat <- matrix(test[], nrow=test@nrows, byrow=TRUE)
image(mat)
persp3d(z = mat, clab = "m")
persp3d(z = mat, col = rainbow(10),border = "black")
persp3d(z = mat, facets = FALSE, curtain = TRUE)

- userSB

---

如果您在某个地方（例如Dropbox）在线发布并提供链接，那么您获得帮助的可能性要大得多。如果您的人口密度数据在单独的文件中，请一并提供。您是想创建一个三维地图（类似于Google Earth），还是一个基于z值着色的区域地图？ - jlhoward

这里有一个链接：https://www.dropbox.com/sh/r35joqzrq5jnhfp/v_jq90etOM 密度的信息在其中。最终我尝试创建一个地图，多边形根据值Z进行着色 :-) 谢谢 - userSB

Answer 1

这符合您的意愿吗？

library(ggplot2)
library(rgdal)           # for readOGR(...) and spTransform(...)
library(RColorBrewer)    # for brewer.pal(...)

setwd("<directory with shapefile>")
map <- readOGR(dsn=".",layer="Ortsteile")
map <- spTransform(map,CRS=CRS("+init=epsg:4839"))
map.data <- data.frame(id=rownames(map@data), map@data)
map.df   <- fortify(map)
map.df   <- merge(map.df,map.data,by="id")
ggplot(map.df, aes(x=long, y=lat, group=group))+
  geom_polygon(aes(fill=Einwohnerd))+
  geom_path(colour="grey")+
  scale_fill_gradientn(colours=rev(brewer.pal(10,"Spectral")))+
  theme(axis.text=element_blank())+
  labs(title="Berlin Ortsteile", x="", y="")+
  coord_fixed()

解释

这是一个非常基本的使用R中的ggplot制作地图的例子，它提供了一个很好的问题。

可以使用readOGR(...)将shapefile读入R中，生成SpatialDataFrame对象。后者基本上有两个部分：一个polygons部分包含多边形边界的坐标，另一个data部分包含来自shapefile属性表中的信息。可以分别引用这些部分，作为map@polygons和map@data。

以上代码读取shapefile并将坐标转换为epsg:4839。然后我们将多边形id（存储在行名中）预置到map@data中的其他信息中，创建map.data。然后我们使用ggplot中的fortify(...)函数将多边形转换为适合绘图的数据框（map.df）。该数据框具有一个id列，该列对应于map.data中的id列。然后我们根据id列将属性信息（map.data）合并到map.df中。

ggplot调用创建地图层并呈现地图，如下所示：

ggplot:       set the default dataset to map.df; identify x- and y-axis columns
geom_polygon: identify column for fill (color of polygon)
geom_path:    polygon boundaries
theme:        turn off axis text
labs:         title, turn off x- and y-axis labels
coord_fixed:  ensures that the map is not distorted

关于scale_fill_gradientn(...)的说明：该函数通过插值colours=参数中提供的调色板来为填充值分配颜色。这里我们使用www.colorbrewer.org的Spectral调色板。不幸的是，该调色板的颜色相反（蓝色-红色），因此我们使用rev(...)来颠倒颜色顺序（高=红色，低=蓝色）。如果你更喜欢matlab中常见的饱和度较高的颜色，请使用library(colorRamps)并将调用scale_fill_gradientn(...)替换为：

  scale_fill_gradientn(colours=matlab.like(10))+

- jlhoward