I'm having trouble adding a shapefile to my ggmap due to differing geographic units

Question

I am trying to add an ESRI shapefile (.shp) to a ggmap plot of the state of North Carolina, which has the following code:

x <- get_map(location="North Carolina", zoom=6, maptype="terrain")

ggmap(x) + geom_polygon(data=citylim83.df, aes(x=long, y=lat), fill="red", alpha=0.2)

The shapefile I've loaded and fortified into citylim83.df . Here is the code used to load the shapefile into ggplot:

    epsgs <- make_EPSG()
    citylim <- readOGR(dsn=".", layer="MunicipalBoundaries_polys")`

The units of the projection of MunicipalBoundaries, after doing a search in EPSG, are ft-US for the state-plane system. Even though this .shp has a geographic coordinate system of NAD83, I want to also project it to NAD83 to get rid of the state plane system (the EPSG code for NAD83 (UTM-17N) is 26917):

    citylim83 <- spTransform(citylim, CRS("+init=epsg:26917"))
    summary(citylim83)

    Object of class SpatialPolygonsDataFrame
    Is projected: TRUE 
    [+init=epsg:26917 +proj=utm +zone=17 +ellps=GRS80 +datum=NAD83 +units=m

    citylim83.df <- fortify(citylim83)

This data frame was then used in the ggmap code shown above.

Even though it is now definitely projected into NAD83, it still will not show up on the base ggmap. What is the base projection of the get_map object I've imported? Is there a command to find this out so I can match my map with the shapefile I want displayed on top of it? Do I have to "unproject" my citylim object? FYI the shapefile is the city limit boundary of each municipality in North Carolina, if it wasn't clear. Any help would be much appreciated on this, as I'm very new to the ggplot2/ggmap community.

Answer 1

Here You have an answer:

library(ggmap)
library(maptools)
shapefile <- readShapeSpatial('MunicipalBoundaries_polys.shp', proj4string = CRS("+proj=lcc +lat_1=34.33333333333334 +lat_2=36.16666666666666 +lat_0=33.75 +lon_0=-79 +x_0=609601.2199999997 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs"))
shp <- spTransform(shapefile, CRS("+proj=longlat +datum=WGS84"))
data <- fortify(shp)

nc <- get_map("North Carolina", zoom = 6, maptype = 'terrain')
ncmap <- ggmap(nc,  extent = "device")
ncmap +
  geom_polygon(aes(x = long, y = lat, group = group), data = data,
               colour = 'grey', fill = 'black', alpha = .4, size = .1)

The coordinates needs to be transform to longlat projection and it's working.