[GIS] ggmap: plot polygon from shapefile

ggmaprshapefile

Using ggmap, I want to include municipality borders (polygon) from a shapefile on a map containing some location points. This script does everything except plotting the polygon:

library(rgdal)
library(ggmap)

# Get shapefile with Drammen municipality borders
tmpzip<-tempfile()
tmpdir<-tempfile()
dir.create(tmpdir)
download.file("http://www.kartverket.no/Documents/Kart/N50-N5000%20Kartdata/33_N5000_shape.zip",tmpzip)
unzip(tmpzip, exdir=tmpdir)
kommune <- readOGR(dsn=tmpdir, layer="NO_AdminOmrader_pol")
kommune<-kommune[kommune$NAVN=="Drammen",]
kommune<-spTransform(kommune, CRS("+init=epsg:4326"))

# Get location point data 
subscr<-data.frame(lon=c(10.1237,10.2161,10.2993),lat=c(59.7567,59.7527,59.6863), pop=c(58,12,150))
coordinates(subscr)<-~lon+lat
proj4string(subscr)<-CRS("+init=epsg:4326")

lon <- c(10.0937,10.3293)
lat <- c(59.7916,59.6563)
map <- get_map(location = c(lon[1], lat[2], lon[2], lat[1]),
               maptype = "roadmap", source = "osm", zoom = 11)
p <- ggmap(map) +
  geom_point(data = as.data.frame(subscr), aes(x = lon, y = lat, size=pop),
             colour = "darkgreen") +
  theme_bw()
print(p)

How can I plot the polygon from the shapefile? I have tried substituting the second last line with the following:

p <- ggmap(map) +
  geom_point(data = as.data.frame(subscr), aes(x = lon, y = lat, size=pop),
             colour = "darkgreen") +
  geom_polygon(data = as.data.frame(kommune)) +
  theme_bw()

But then I get the following error:

Error: Aesthetics must be either length 1 or the same as the data (1): x, y

Best Answer

as.data.frame() does not work for SpatialPolgons in geom_polygon, because the geometry gets lost. You have to use ggplot2::fortify (may be deprecated in the future, see ?fortify). The recommended way is now to use broom::tidy:

R> library("broom")
R> head(tidy(kommune))
Regions defined for each Polygons
   long   lat order  hole piece group  id
1 10.29 59.72     1 FALSE     1 153.1 153
2 10.32 59.70     2 FALSE     1 153.1 153
3 10.32 59.69     3 FALSE     1 153.1 153
4 10.31 59.68     4 FALSE     1 153.1 153
5 10.30 59.67     5 FALSE     1 153.1 153
6 10.28 59.67     6 FALSE     1 153.1 153

But another problem arises with your example. Since the polygon is larger than the map extent, ggmap does not correctly clip the polygon. ggmap sets the limits on the scale, this will throw away all data that's not inside these limits.

Here is a modified version of your code:

p <- ggmap(map, extent = "normal", maprange = FALSE) +
     geom_point(data = as.data.frame(subscr),
                aes(x = lon, y = lat, size=pop),
                colour = "darkgreen") +
     geom_polygon(data = fortify(kommune),
                  aes(long, lat, group = group),
                  fill = "orange", colour = "red", alpha = 0.2) +
     theme_bw() +
     coord_map(projection="mercator",
               xlim=c(attr(map, "bb")$ll.lon, attr(map, "bb")$ur.lon),
               ylim=c(attr(map, "bb")$ll.lat, attr(map, "bb")$ur.lat))

print(p)

ggmap

Related Question