R overlay geom_polygon with ggmap object, spatial file conversion

Question

I've seen examples of this in several places including in the ggmap article in the Rjournal (https://journal.r-project.org/archive/2013-1/kahle-wickham.pdf) and see this for another walkthrough - https://markhneedham.com/blog/2014/11/17/r-ggmap-overlay-shapefile-with-filled-polygon-of-regions/

The problem I'm having is simply implementing this. It seems straight forward but I am missing something.

I am using a shape file of Wisconsin Counties from the Wisconsin Dept. of Natural Resources (free) https://data-wi-dnr.opendata.arcgis.com/datasets/8b8a0896378449538cf1138a969afbc6_3?geometry=-110.743%2C42.025%2C-68.93%2C47.48

Here's the code:

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

I can plot the shapefile just fine using plot(shpfile). Next I convert this to a format suitable for plotting in ggplot. Many examples use "fortify" but that seems to have been replaced by "tidy" which is part of the "broom" package. FWIW, I've tried it with fortify and get the same result.

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

tidydta <- tidy(shpfile, group=group)

Now I can successfully plot the shapefile in ggplot as a polygon

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

Next I retrieve a background map using ggmap.

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

The problem is that I can't combine these. I'm assuming there must be something wrong with the tidy conversion or I'm missing a step. I do get an error:

in min(x): no non-missing arguments to min; returning Inf

which I think happens because I have a zero length vector somewhere.

Here's the command:

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

I have successfully added geo-coded points to the map using geom_point but I'm stuck with the polygon.

Can anyone tell me what I'm doing wrong?

Z.Lin Z.Lin · Accepted Answer · 2019-03-08T10:22:24

The coordinate system used by your shapefile isn't lat-lon. You should transform it before converting it into a dataframe for ggplot. The following works:

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)

For future reference, do check your coordinate values by printing the dataframe to console / plotting with visible x/y axis labels. If the numbers are wildly different from that of your background map's bounding box (e.g. 7e+05 vs. 47), you probably need to make some transformations. 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")
)

R overlay geom_polygon with ggmap object, spatial file conversion

1 Answers