Я пытался использовать данные переписи, которые вы имеете в своем сценарии. Однако R Studio почему-то продолжал падать, когда я применил read.shp() к данным многоугольника. Поэтому я решил использовать пример со страницы справки read.shp(), которая также является данными переписи. Я надеюсь, ты не против. Потребовалось некоторое время, чтобы понять, как нарисовать карту с классом shp. Позвольте мне объяснить, что я прошел шаг за шагом.
Эта часть со страницы справки. Я в основном получаю шейп-файл и импортирую его как объект 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 ...
Вот пример для Аляски. Как вы видите, в * 1016 есть некоторые числа * Эти числа указывают, где начинаются новые данные многоугольника. Алакса имеет много небольших островов. Следовательно, они должны были указать разные полигоны в данных с этой информацией. Мы вернемся к этому позже, когда создадим фреймы данных.
#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, чтобы мы могли указывать все полигоны с уникальными идентификаторами. Нам нужно создать новую групповую переменную, которая содержит уникальное значение идентификатора для каждого многоугольника. Я использовал findInterval(), который использует индексы для создания групповой переменной. Одна сложность заключалась в том, что нам нужно использовать left.open = TRUE в findInterval() для создания групповой переменной. (Это дало мне некоторое время, чтобы выяснить, что происходит.) Эта часть 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 dr aws забавные длинные строки, которые вы можете увидеть даже в своем примере. Нам нужно преобразовать эти положительные числа в отрицательные, чтобы 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))
