使用带有 ggplot2 的 read.shp 绘制形状文件

Question

我想绘制一个使用 fastshp 包中的 read.shp 加载的形状文件。但是，read.shp 函数返回列表列表而不是 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() 帮助页面中的示例，这也是人口普查数据。希望你不要介意。花了一些时间才弄清楚如何使用 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))