Add shapefile to ggmap

Question

Somewhat similar to the issue described here I am having troubles aligning a shapefile and ggmap object.

1. My shapefile consists of local area boundaries in Victoria, Australia, and I am trying to overlay them on top of a google map of the state (Victoria).

   The source shapefile has the following PROJ4 string (extracted from the prj file)

   [+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs]
   

   which corresponds to EPSG:4283.

   Here's the summary of my shapefile object sp:

   > summary(sp)
   Object of class SpatialPolygonDataFrame
   Coordinates:
           min        max
   x  96.81677 159.109219
   y -43.74051  -9.142176
   Is projected: FALSE
   proj4string:
   [+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs]
   

2. I convert the shapefile to match google's pseudo Mercator projection (at least that's what I think I'm doing)

   sp <- spTransform(sp, CRS("+proj=longlat +init=epsg:3857"))
   

   and turn sp into a fortified dataframe df.sp.

3. I then use

   map <- get_map("Victoria", zoom = 7, maptype = "terrain", source = "google")
   

   to get a google terrain map of Victoria, and (gg)plot them

   ggmap(map) + geom_polygon(data = df.sp, aes(x = long, y = lat, group = group)) + coord_equal() + theme_map()
   

It's clear from the resulting plot that shapefile coordinates and googlemap coordinates don't overlap. Am I doing something wrong with the coordinate transformation? How do I properly match shapefile and googlemap coordinates? I would appreciate any help/insight into this matter.

Answer 1

This is my first post, so forgive me if the styling isn't fantastic. I found for whatever reason, the default ESPG:3857 in the rgdal package wasn't transforming the projection correctly in the spTransform function.

• Using gene's wonderful answer in GIS StackExchange, I was directed to the Prj2EPSG website and the Spatial Reference website. I double checked the ABS's proj4 string (I assume that's where this data has come from) and it's as you've said (which is EPSG:4283).

• When entering EPSG:3857 into the Spatial reference I got this string (using the Proj4 link) +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs.

• I replaced +proj=merc with +proj=longlat to project the longitude and latitude co-ordinates.

• Using that string with merc replaced seems to give me better results than using just +init=epsg:3857, which indicated to me there's some difference in the way the rgdal package handles it.

I'd love to hear of a way to improve on this further, but I'm hoping maybe this can help you get better results in the meantime.

---

EDIT: I've also found that when plotting using ggplot2 adding coord_map() seems to fix the shape files' alignment when up close, but leave them a little off when zoomed out further. At zoom scales 7 and above I seem to get perfect alignment, but 6 and lower and I don't. It seems that for whatever reason the different zoom levels change the projection slightly.

So try

ggmap(map) + 
  geom_polygon(data = df.sp, aes(x = long, y = lat, group = group)) +
  coord_equal() + 
  theme_map() + 
  coord_map()