R наложение geom_polygon с объектом ggmap, пространственное преобразование файлов - PullRequest
Новая
       35

R наложение geom_polygon с объектом ggmap, пространственное преобразование файлов

2 голосов
/ 08 марта 2019

Я видел примеры этого в нескольких местах, в том числе в статье ggmap в журнале (https://journal.r -project.org / archive / 2013-1 / kahle-wickham.pdf ) и посмотрите это для другого прохождения - https://markhneedham.com/blog/2014/11/17/r-ggmap-overlay-shapefile-with-filled-polygon-of-regions/

Проблема, с которой я столкнулся, заключается в простой реализации этого. Это кажется прямым, но я что-то упускаю.

Я использую файл формы округов Висконсин из Висконсинского департамента природных ресурсов (бесплатно) https://data -wi-dnr.opendata.arcgis.com / наборов данных / 8b8a0896378449538cf1138a969afbc6_3? Геометрия = -110,743% 2C42.025% 2C-68,93% 2C47.48

Вот код: 

library(rgdal)
shpfile <- readOGR(dsn = "[file path to the shapefile directory]", 
                   stringsAsFactors = FALSE )

Я могу нарисовать шейп-файл, используя plot(shpfile). Далее я конвертирую это в формат, подходящий для печати в ggplot. Во многих примерах используется «fortify», но, похоже, его заменили на «tidy», который является частью пакета «broom». FWIW, я попробовал это с Fortify и получить тот же результат. 

library(broom)
library(ggplot2)
library(ggmap)

tidydta <- tidy(shpfile, group=group)

Теперь я могу успешно построить шейп-файл в ggplot в виде многоугольника 

ggplot() +
  geom_polygon(data=tidydta, 
               mapping=aes(y=lat , x=long, group=group), 
               color="dark red", alpha=.2) +
  theme_void()

Далее я получаю фоновую карту, используя ggmap. 

wisc <- get_map(center = c(lon= -89.75, lat=44.75), zoom=7, maptype="toner")

Проблема в том, что я не могу объединить их. Я предполагаю, что должно быть что-то не так с аккуратным преобразованием, или я пропускаю шаг. Я получаю сообщение об ошибке:

в мин (х): нет не пропущенных аргументов мин; возвращая Inf

что, как мне кажется, происходит потому, что у меня где-то есть вектор нулевой длины.

Вот команда: 

ggmap(wisc) +
  geom_polygon(aes(x=long, y=lat, group=group), 
               data=tidydta, 
               color="dark red", alpha=.2, size=.2)

Я успешно добавил геокодированные точки на карту, используя geom_point, но я застрял с многоугольником.

Может кто-нибудь сказать мне, что я делаю не так?

1 Ответ

2 голосов
/ 08 марта 2019

Система координат, используемая вашим шейп-файлом, не широка. Вы должны преобразовать его перед преобразованием в фрейм данных для ggplot. Следующие работы: 

shpfile <- spTransform(shpfile, "+init=epsg:4326") # transform coordinates
tidydta2 <- tidy(shpfile, group=group) 

wisc <- get_map(location = c(lon= -89.75, lat=44.75), zoom=7)

ggmap(wisc) +
  geom_polygon(aes(x=long, y=lat, group=group), 
               data=tidydta2, 
               color="dark red", alpha=.2, size=.2)

plot

Для справки в будущем, проверьте значения координат, распечатав кадр данных на консоли / построив график с видимыми метками оси x / y. Если числа сильно отличаются от чисел в ограничительной рамке вашей фоновой карты (например, 7e + 05 против 47), вам, вероятно, необходимо выполнить некоторые преобразования. E.g.:

> head(tidydta) # coordinate values of dataframe created from original shapefile
# A tibble: 6 x 7
     long     lat order hole  piece group id   
    <dbl>   <dbl> <int> <lgl> <chr> <chr> <chr>
1 699813. 246227.     1 FALSE 1     0.1   0    
2 699794. 246082.     2 FALSE 1     0.1   0    
3 699790. 246007.     3 FALSE 1     0.1   0    
4 699776. 246001.     4 FALSE 1     0.1   0    
5 699764. 245943.     5 FALSE 1     0.1   0    
6 699760. 245939.     6 FALSE 1     0.1   0    

> head(tidydta2) # coordinate values of dataframe created from transformed shapefile
# A tibble: 6 x 7
   long   lat order hole  piece group id   
  <dbl> <dbl> <int> <lgl> <chr> <chr> <chr>
1 -87.8  42.7     1 FALSE 1     0.1   0    
2 -87.8  42.7     2 FALSE 1     0.1   0    
3 -87.8  42.7     3 FALSE 1     0.1   0    
4 -87.8  42.7     4 FALSE 1     0.1   0    
5 -87.8  42.7     5 FALSE 1     0.1   0    
6 -87.8  42.7     6 FALSE 1     0.1   0 

> attr(wisc, "bb") # bounding box of background map
# A tibble: 1 x 4
  ll.lat ll.lon ur.lat ur.lon
   <dbl>  <dbl>  <dbl>  <dbl>
1   42.2  -93.3   47.2  -86.2

# look at the axis values; don't use theme_void until you've checked them
cowplot::plot_grid(
  ggplot() +
    geom_polygon(aes(x=long, y=lat, group=group), 
                 data=tidydta, 
                 color="dark red", alpha=.2, size=.2),
  ggplot() +
    geom_polygon(aes(x=long, y=lat, group=group), 
                 data=tidydta2, 
                 color="dark red", alpha=.2, size=.2),
  labels = c("Original", "Transformed")
)

comparison

...