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

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

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

Я использую шейп-файл округов Висконсин из Департамента природных ресурсов Висконсина (бесплатно) https://data-wi-dnr.opendata.arcgis.com/datasets/8b8a0896378449538cf1138a969afbc6_3%geometry=-110.7432C42.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")

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

in min (x): для min нет непустых аргументов; возвращение Inf

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

Вот команда:

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

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

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


r ggplot2 shapefile ggmap spatial
person Will Hauser    schedule 07.03.2019    source источник

Ответы (1)


arrow_upward
3
arrow_downward

Система координат, используемая вашим шейп-файлом, не является широтой и долготой. Вы должны преобразовать его перед преобразованием во фрейм данных для 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)

сюжет

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

> 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")
)

сравнение

person Z.Lin    schedule 08.03.2019
comment
Бесконечно благодарен! Я видел, как sptransform используется здесь и там, но не видел краткого объяснения того, как определить, когда это необходимо. Это очень полезно. - person Will Hauser; 08.03.2019