绘制使用 read.shp 和 ggplot2 加载的形状文件

Question

我想绘制一个使用 fastshp 包中的 read.shp 加载的形状文件。但是，read.shp 函数是 returns 列表而不是 data.frame 列表。我不确定我需要提取列表的哪一部分才能获得格式正确的 data.frame 对象。 stack overflow 上已经有人问过这个确切的问题，但是，该解决方案似乎不再有效（解决方案来自 > 7 年前）。非常感谢任何帮助。

remotes::install_github("s-u/fastshp") #fastshp not on CRAN
library(ggplot2);library(fastshp)

temp <- tempfile()
temp2 <- tempfile()
download.file("https://www2.census.gov/geo/tiger/TIGER2017/COUNTY/tl_2017_us_county.zip",temp)
unzip(zipfile = temp, exdir = temp2)
shp <- list.files(temp2, pattern = ".shp$",full.names=TRUE) %>% read.shp(.)

shp 是包含大量信息的列表列表。我尝试了之前发布的 SO 中的以下解决方案，但无济于事：

shp.list <- sapply(shp, FUN = function(x) Polygon(cbind(lon = x$x, lat = x$y))) #throws an error here cbind(lon = x$x, lat = x$y) returns NULL
shp.poly <- Polygons(shp.list, "area")
shp.df <- fortify(shp.poly, region = "area")

我还尝试了以下方法：

shp.list <- sapply(shp, FUN = function(x) do.call(cbind, x[c("id","x","y")])) #returns NULL value here...
shp.df <- as.data.frame(do.call(rbind, shp.list))

已更新：仍然没有运气但更接近：

file_shp<-list.files(temp2, pattern = ".shp$",full.names=TRUE) %>%
  read.shp(., format = c("table"))

ggplot() + 
geom_polygon(data = file_shp, aes(x = x, y = y, group = part), 
             colour = "black", fill = NA)

投影似乎关闭了。我不确定如何命令数据正确映射，也不确定如何读取 CRS 数据。尝试了以下无济于事：

file_prj<-list.files(temp2, pattern = ".prj$",full.names=TRUE) %>%
  proj4string(.)

Answer 1

我尝试使用您脚本中的人口普查数据。但是，当我将 read.shp() 应用于多边形数据时，R Studio 不知何故一直崩溃。因此，我决定使用 read.shp() 帮助页面中的示例，该示例也是人口普查数据。希望你不会介意。花了一些时间弄清楚如何用 class shp 绘制地图。让我一步一步解释一下我经历了什么。

这部分来自帮助页面。我基本上是在获取 shapefile 并将其作为 shp 对象导入。

# Census 2010 TIGER/Line(TM) state shapefile
library(fastshp)
fn <- system.file("shp", "tl_2010_us_state10.shp.xz", package="fastshp")
s <- read.shp(xzfile(fn, "rb"))

让我们检查一下这个对象，s 是什么样的。它包含 52 个列表。在每个列表中，有六个向量。 ID 是一个唯一的整数来表示一个状态。 x 是经度，y 是纬度。最糟糕的部分是 parts。在下面的这个例子中，只有一个数字，这意味着在这种状态下只有一个多边形。但是其他一些列表（状态）有多个数字。这些数字基本上是指示数据中新多边形开始位置的索引。

#> str(s)
#List of 52
# $ :List of 6
#  ..$ id   : int 1
#  ..$ type : int 5
#  ..$ box  : num [1:4] -111 41 -104 45
#  ..$ parts: int 0
#  ..$ x    : num [1:9145] -109 -109 -109 -109 -109 ...
#  ..$ y    : num [1:9145] 45 45 45 45 45 ...

这是阿拉斯加的。如您所见，parts 中有一些数字，这些数字表示新多边形数据的开始位置。阿拉克萨有许多小岛。因此，他们需要使用此信息指示数据中的不同多边形。稍后我们将在创建数据框时回到这一点。

#List of 6
# $ id   : int 18
# $ type : int 5
# $ box  : num [1:4] -179.2 51.2 179.9 71.4
# $ parts: int [1:50] 0 52 88 127 175 207 244 306 341 375 ...
# $ x    : num [1:14033] 177 177 177 177 177 ...
# $ y    : num [1:14033] 52.1 52.1 52.1 52.1 52.1 ...

我们需要的是以下内容。对于每个列表，我们需要提取经度（即 x）、纬度（即 y）和 id，以便为一个州创建数据名。此外，我们需要使用 parts 以便我们可以用唯一 ID 表示所有多边形。我们需要创建一个新的组变量，其中包含每个多边形的唯一 ID 值。我使用 findInterval() 它采用索引来创建组变量。一个棘手的部分是我们需要在 findInterval() 中使用 left.open = TRUE 来创建一个组变量。（这让我很难弄清楚发生了什么。）这个 map_dfr() 部分处理我刚才描述的工作。

library(tidyverse)

map_dfr(.x = s,
        .f = function(mylist){

                temp <- data.frame(id = mylist$id,
                                   lon = mylist$x,
                                   lat = mylist$y)
                ind <- mylist$parts

                out <- mutate(temp,
                              subgroup = findInterval(x = 1:n(), vec = ind, left.open = TRUE),
                              group = paste(id, subgroup, sep = "_"))
                return(out)

                }) -> test

一旦我们有了 test，我们就有了另一份工作。阿拉斯加的一些经度点保持正数（例如，179.85）。只要我们有这样的数字，ggplot2 就会画出有趣的长线，即使在您的示例中也可以看到。我们需要的是将这些正数转换为负数，以便ggplot2可以绘制出合适的地图。

mutate(test,
       lon = if_else(lon > 0, lon * -1, lon)) -> out

此时，out 看起来像这样。

  id       lon      lat subgroup group
1  1 -108.6213 45.00028        1   1_1
2  1 -108.6197 45.00028        1   1_1
3  1 -108.6150 45.00031        1   1_1
4  1 -108.6134 45.00032        1   1_1
5  1 -108.6133 45.00032        1   1_1
6  1 -108.6130 45.00032        1   1_1

现在我们准备绘制地图。

ggplot() +
geom_polygon(data = out, aes(x = lon, y = lat, group = group))

绘制使用 read.shp 和 ggplot2 加载的形状文件

Plotting shape files loaded using read.shp with ggplot2

maps

r

spatial

ggplot2