Plotting OpenStreetMap with ggmap

4
votes

I'm trying to get districts of Warsaw and draw them on google map. Using this code, where 2536107 is relation code for OpenStreetMap single Warsaw district, gives me almost what I want but with a few bugs. There is general outline but also lines between points which shouldn't be connected. What am I doing wrong? 

map <- get_googlemap('warsaw', zoom =10) 
warszawa <- get_osm(relation(2536107), full = T)
warszawa.sp <- as_sp(warszawa, what='lines')
warsawfort <- fortify(warszawa.sp)

mapa_polski <- ggmap(map, extent='device', legend="bottomleft") 
warsawfort2 <- geom_polygon(aes(x = long, y = lat), 
               data = warsawfort, fill="blue", colour="black", 
               alpha=0.0, size = 0.3)

base <- mapa_polski + warsawfort2 
base

Edit: I figured it must be somehow connected with order of plotting every point/line but I have no idea how to fix this.

2
rggplot2openstreetmapggmap

2 Answers

9
votes

There is a way to generate your map without using external packages: don't use osmar...

This link, to the excellent Mapzen website, provides a set of shapefiles of administrative areas in Poland. If you download and unzip it, you will see a shapfile set called warsaw.osm-admin.*. This is a polygon shapefile of all the districts in Poland, conveniantly indexed by osm_id(!!). The code below assumes you have downloaded the file and unzipped it into the "directory with your shapefiles".

library(ggmap)
library(ggplot2)
library(rgdal)

setwd(" <directory with your shapefiles> ")
pol    <- readOGR(dsn=".",layer="warsaw.osm-admin")
spp    <- pol[pol$osm_id==-2536107,]
wgs.84 <- "+proj=longlat +datum=WGS84"
spp    <- spTransform(spp,CRS(wgs.84))

map    <- get_googlemap('warsaw', zoom =10) 
spp.df <- fortify(spp)

ggmap(map, extent='device', legend="bottomleft") +
  geom_polygon(data = spp.df, aes(x = long, y=lat, group=group), 
               fill="blue", alpha=0.2) +
  geom_path(data=spp.df, aes(x=long, y=lat, group=group), 
            color="gray50", size=0.3)

Two nuances: (1) The osm IDs are stored as negative numbers, so you have to use, e.g.,

spp    <- pol[pol$osm_id==-2536107,]

to extract the relevant district, and (2) the shapefile is not projected in WGS84 (long/lat). So we have to reproject it using:

spp    <- spTransform(spp,CRS(wgs.84))

The reason osmar doesn't work is that the paths are in the wrong order. Your warszawa.sp is a SpatialLinesDataframe, made up of a set of paths (12 in your case), each of which is made up of a set of line segments. When you use fortify(...) on this, ggplot tries to combine them into a single sequence of points. But since the paths are not in convex order, ggplot tries, for example, to connect a path that ends in the northeast, to a path the begins in the southwest. This is why you're getting all the extra lines. You can see this by coloring the segments:

xx=coordinates(warszawa.sp)
colors=rainbow(11)
plot(t(bbox(warszawa.sp)))
lapply(1:11,function(i)lines(xx[[i]][[1]],col=colors[i],lwd=2))

The colors are in "rainbow" order (red, orange, yellow, green, etc.). Clearly, the lines are not in that order.

EDIT Response to @ako's comment.

There is a way to "fix" the SpatialLines object, but it's not trivial. The function gPolygonize(...) in the rgeos package will take a list of SpatialLines and convert to a SpatialPolygons object, which can be used in ggplot with fortify(...). One huge problem (which I don't understand, frankly), is that OP's warszaw.sp object has 12 lines, two of which seem to be duplicates - this causes gPolygonize(...) to fail. So if you create a SpatialLines list with just the first 11 paths, you can convert warszawa.sp to a polygon. This is not general however, as I can't predict how or if it would work with other SpatialLines objects converted from osm. Here's the code, which leads to the same map as above.

library(rgeos)
coords <- coordinates(warszawa.sp)
sll <- lapply(coords[1:11],function(x) SpatialLines(list(Lines(list(Line(x[[1]])),ID=1))))
spp <- gPolygonize(sll)
spp.df <- fortify(spp)
ggmap(map, extent='device', legend="bottomleft") +
  geom_polygon(data = spp.df, aes(x = long, y=lat, group=group), 
               fill="blue", alpha=0.2) +
  geom_path(data=spp.df, aes(x=long, y=lat, group=group), 
            color="gray50", size=0.3)
2
votes

I am not sure this is a general hangup--I can reproduce your example and see the issue. My first thought was that you didn't supply group=id which are typically used for polygons with many lines, but you have lines, so that should not be needed.

The only way I could get it to display properly was by changing your lines into a polygon off script. Qgis' line to polygon didn't get this "right", getting a large donut hole, so I used ArcMap, which produced a full polygon. If this is a one off that may work for your workflow. Odds are it is not. In that case, perhaps RGDAL can transform lines to polygons, assuming that is indeed a general problem.

Upon reading the polygon shapefile and fortifying that, your code ran without problems. enter image description here