Choropleth map in ggplot with polygons that have holes

Question

I'm trying to draw a choropleth map of Germany showing poverty rate by state (inspired by this question).

The problem is that some of the states (Berlin, for example) are completely surrounded by other states (Brandenburg), and I'm having trouble getting ggplot to recognize the "hole" in Brandenburg.

The data for this example is here.

library(rgdal)
library(ggplot2)
library(RColorBrewer)

map <- readOGR(dsn=".", layer="germany3")
pov <- read.csv("gerpoverty.csv")

mrg.df <- data.frame(id=rownames(map@data),ID_1=map@data$ID_1)
mrg.df <- merge(mrg.df,pov, by="ID_1")
map.df <- fortify(map)
map.df <- merge(map.df,mrg.df[,c("id","poverty")], by="id")
ggplot(map.df, aes(x=long, y=lat, group=group)) +
  geom_polygon(aes(fill=poverty))+
  geom_path(colour="grey50")+
  scale_fill_gradientn(colours=brewer.pal(5,"OrRd"))+
  labs(x="",y="")+ theme_bw()+
  coord_fixed()

Notice how the colors for Berlin and Brandenburg (in the northeast) are identical. They shouldn't be - Berlin's poverty rate is much lower than Brandenburg. It appears that ggplot is rendering the Berlin polygon and then rendering the Brandenburg polygon over it, without the hole.

If I change the call to geom_polygon(...) as suggested here, I can fix the Berlin/Brandenburg problem, but now the three northernmost states are rendered incorrectly.

ggplot(map.df, aes(x=long, y=lat, group=group)) +
  geom_polygon(aes(group=poverty, fill=poverty))+
  geom_path(colour="grey50")+
  scale_fill_gradientn(colours=brewer.pal(5,"OrRd"))+
  labs(x="",y="")+ theme_bw()+
  coord_fixed()

What am I doing wrong??

Reverberation

Did you try to use map<-fortify(map) for your map? docs.ggplot2.org/0.9.3.1/fortify.sp.html – Occasionally

See line 8 of the code: map.df <- fortify(map). Or did you mean something else? – Reverberation

There is a discussion and example of a workaround for this problem at github.com/hadley/ggplot2/wiki/plotting-polygon-shapefiles – Theroid

@Theroid - Thanks for this. Do you know if this problem is going to be fixed? Other packages don't fail this way (see the answer below). Also, I do not think this will work if there is "just" a hole (e.g. a lake) which should render as transparent (background shows through). This response has an example of this. – Reverberation

@Reverberation - I doubt this will be fixed anytime soon, you'll have to workaround as I've illustrated. – Theroid

Answer 1

You can plot the island polygons in a separate layer, following the example on the ggplot2 wiki. I've modified your merging steps to make this easier:

mrg.df <- data.frame(id=rownames(map@data),ID_1=map@data$ID_1)
mrg.df <- merge(mrg.df,pov, by="ID_1")
map.df <- fortify(map)
map.df <- merge(map.df,mrg.df, by="id")

ggplot(map.df, aes(x=long, y=lat, group=group)) +
    geom_polygon(aes(fill=poverty), color = "grey50", data =subset(map.df, !Id1 %in% c("Berlin", "Bremen")))+
    geom_polygon(aes(fill=poverty), color = "grey50", data =subset(map.df, Id1 %in%  c("Berlin", "Bremen")))+
    scale_fill_gradientn(colours=brewer.pal(5,"OrRd"))+
    labs(x="",y="")+ theme_bw()+
    coord_fixed()

As an unsolicited act of evangelism, I encourage you to consider something like

library(ggmap)
qmap("germany", zoom = 6) +
    geom_polygon(aes(x=long, y=lat, group=group, fill=poverty),
                 color = "grey50", alpha = .7,
                 data =subset(map.df, !Id1 %in% c("Berlin", "Bremen")))+
    geom_polygon(aes(x=long, y=lat, group=group, fill=poverty),
                 color = "grey50", alpha= .7,
                 data =subset(map.df, Id1 %in%  c("Berlin", "Bremen")))+
    scale_fill_gradientn(colours=brewer.pal(5,"OrRd"))

to provide context and familiar reference points.

Theroid

Answer 2

This is just an expansion on @Ista's answer, which does not require that one knows which states (Berlin, Bremen) need to be rendered last.

This approach takes advantage of the fact that fortify(...) generates a column, hole which identifies whether a group of coordinates are a hole. So this renders all regions (id's) with any holes before (e.g. underneath) the regions without holes.

Many thanks to @Ista, without whose answer I could not have come up with this (believe me, I spent many hours trying...)

ggplot(map.df, aes(x=long, y=lat, group=group)) +
  geom_polygon(data=map.df[map.df$id %in% map.df[map.df$hole,]$id,],aes(fill=poverty))+
  geom_polygon(data=map.df[!map.df$id %in% map.df[map.df$hole,]$id,],aes(fill=poverty))+
  geom_path(colour="grey50")+
  scale_fill_gradientn(colours=brewer.pal(5,"OrRd"))+
  labs(x="",y="")+ theme_bw()+
  coord_fixed()

Reverberation

@jhoward - I actually tried and failed to do this for my original answer. I had thought that 'geom_polygon(data=map.df[!map.df$hole,],aes(fill=poverty))' should work and was frustrated that it didn't. Nice to see your solution here. – Theroid

@jhoward Thanks a million for this brilliant solution! – Moneylender

Answer 3

Just to add another small improvement to @Ista's and @jhoward's answers (thanks a lot for your help!).

The modification of @jhoward could be easily wrapped in a small function like this

gghole <- function(fort){
        poly <- fort[fort$id %in% fort[fort$hole,]$id,]
        hole <- fort[!fort$id %in% fort[fort$hole,]$id,]
        out <- list(poly,hole)
        names(out) <- c('poly','hole')
        return(out)
} 
# input has to be a fortified data.frame

Then, one doesn't need to recall every time how to extract holes info. The code would look like

    ggplot(map.df, aes(x=long, y=lat, group=group)) +
            geom_polygon(data=gghole(map.df)[[1]],aes(fill=poverty),colour="grey50")+
            geom_polygon(data=gghole(map.df)[[2]],aes(fill=poverty),colour="grey50")+
    # (optionally). Call by name
    #         geom_polygon(data=gghole(map.df)$poly,aes(fill=poverty),colour="grey50")+
    #         geom_polygon(data=gghole(map.df)$hole,aes(fill=poverty),colour="grey50")+
            scale_fill_gradientn(colours=brewer.pal(5,"OrRd"))+
            labs(x="",y="")+ theme_bw()+
            coord_fixed()

Moneylender

a shorter version: ggpolyhole <- function(fort, poly = TRUE){ idx <- fort$id %in% fort[fort$hole, ]$id fort[idx - poly, ] } – Soloman

Answer 4

Alternatively you could create that map using rworldmap.

library(rworldmap)
library(RColorBrewer)
library(rgdal)

map <- readOGR(dsn=".", layer="germany3")
pov <- read.csv("gerpoverty.csv")

#join data to the map
sPDF <- joinData2Map(pov,nameMap='map',nameJoinIDMap='VARNAME_1',nameJoinColumnData='Id1')

#default map
#mapPolys(sPDF,nameColumnToPlot='poverty')

colours=brewer.pal(5,"OrRd")
mapParams <- mapPolys( sPDF
                      ,nameColumnToPlot='poverty'
                      ,catMethod="pretty"
                      ,numCats=5
                      ,colourPalette=colours
                      ,addLegend=FALSE )


do.call( addMapLegend, c( mapParams
                          , legendLabels="all"
                          , legendWidth=0.5
                        ))

#to test state names
#text(pov$x,pov$y,labels=pov$Id1)

Armijo

Thanks, but it's not about this specific map. I'm trying to find out if this is a bug in ggplot, or if I'm doing something wrong. – Reverberation