Add shapefile to ggmap

Question

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

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]
```
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.

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.

This doesn't answer your question, but I've got an example of plotting Victorian polygons on a Google Map in the vignette of my googleway package — SymbolixAU

LachlanO LachlanO · Accepted Answer · 2017-12-20T01:21:20

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

Add shapefile to ggmap

1 Answers