Region polygons not displaying in ggplot2 Choropleth map

Question

I'm trying to do a very basic map plot with ggplot2. I can't find why the colored polygons would not display. It seems that my code is no different from what I can find on the many tutorials and the questions already answered on this website. I think it could come from the way I prepare the data (see 100% reproducible example below).

library(maptools)
library(sp)
library(ggplot2)
library(rgeos)
con <- url("http://biogeo.ucdavis.edu/data/gadm2/R/GHA_adm1.RData")
print(load(con))
close(con)

ghaDF<-as.data.frame(gadm)
ghaDF$prod <- c(12, 26, 12,22,0,11,4,5,4,4) # add values for the regions to be colored with
gadm@data$id = rownames(gadm@data) #create an id in the shapefile data
ghaMap <- fortify(gadm, region="id")
colnames(ghaDF[5])<-"id"
ghaMap <- merge(ghaMap, ghaDF)

m0 <- ggplot(data=ghaMap)
m1 <- m0 + geom_polygon(aes(x=long, y=lat, fill = prod, group=group)) 
         + scale_fill_gradient(low = "light green", high = "dark green")
m2 <- m1 + geom_path(aes(x=long, y=lat, group=group),color='gray')  
         + coord_equal()
m2

On the image above (the output of m2), the regions should be colored according to the ghaMap$prod variable. Any suggestions?

(R version 3.0.2 - Platform: x86_64-w64-mingw32/x64 (64-bit))

Answer 1

Your data seems fine. Here's a solution with geom_map (and, colorbrewer colors :-):

devtools::source_gist("https://gist.github.com/hrbrmstr/33baa3a79c5cfef0f6df")
gg <- ggplot()
gg <- gg + geom_map(data=ghaMap, map=ghaMap,
                    aes(map_id=id, group=group,
                        x=long, y=lat, fill=prod))
gg <- gg + scale_fill_gradient(low="#99d8c9", high="#00441b")
gg <- gg + coord_map()
gg <- gg + theme_map() + theme(legend.position="right")
gg

You can omit the devtools… and theme_map() without impacting the solution.

EDIT Using simplified polygons

library(maptools)
library(sp)
library(ggplot2)
library(rgeos)
library(rgdal)

# get the shapefile
download.file("http://biogeo.ucdavis.edu/data/gadm2/shp/GHA_adm.zip", "GHA_adm.zip")
unzip("GHA_adm.zip", exdir="GHA_adm")

# simplify it
setwd("GHA_adm")
system("ogr2ogr -simplify 0.001 simple1.shp GHA_adm1.shp") # simplify (500K -> 80K)
setwd("..")

# read it in
gadm <- readOGR("GHA_adm", "simple1")

# convert it
gadm_map <- fortify(gadm, region="NAME_1")

# make the values we want to fill with
# you can use "ID_1" but folks I've seen generally like using names. works either way
prod <- data.frame(id=gadm$NAME_1, value=c(12, 26, 12, 22, 0, 11, 4, 5, 4, 4), stringsAs=FALSE)

# merge the data together
gadm_map <- merge(gadm_map, prod, all.x=TRUE) # add it right into the fortified data frame

# plot it
devtools::source_gist("https://gist.github.com/hrbrmstr/33baa3a79c5cfef0f6df")
gg <- ggplot()
gg <- gg + geom_map(data=gadm_map, map=gadm_map,
                    aes(map_id=id, group=group,
                        x=long, y=lat, fill=value))
gg <- gg + scale_fill_gradient(low="#99d8c9", high="#00441b")
gg <- gg + coord_map()
gg <- gg + theme_map() + theme(legend.position="right")
gg

(same map as above with slightly less edge detail)

Using lsos (via):

lsos()

##                               Type   Size Rows Columns
## gadm_map                data.frame 275256 5167       9
## gadm      SpatialPolygonsDataFrame 168296   10       9
## gg                              gg  34792    9      NA
## prod                    data.frame   2136   10       3

vs the same from the OP:

##                               Type    Size  Rows Columns
## ghaMap                  data.frame 3689120 31735      18
## gadm      SpatialPolygonsDataFrame  589424    10      11
## gg                              gg   34792     9      NA
## ghaDF                   data.frame    5088    10      12
## con                            url     552     1      NA

Answer 2

Your problem is very simple. Change

colnames(ghaDF[5])<-"id"

to

colnames(ghaDF)[5]<-"id"

and your code produces the map you want.

The first line above extracts column 5 of ghaDF as a data frame with 1 column, sets that column name to "id", and discards it when the call to colnames(...) completes. It does not affect the original data frame at all.

The second line extracts the column names of ghaDF as a character vector, and sets the 5th element of that vector to "id", which in effect changes the name of column 5 of ghaDF.

Having said all this, your workflow is a bit tortured. For one thing you are merging all the columns of gadm@data into ghaMap, which is unnecessary and extremely inefficient. The code below produces the same map:

load(url("http://biogeo.ucdavis.edu/data/gadm2/R/GHA_adm1.RData"))
ghaMap <- fortify(gadm)
ghaDF  <- data.frame(id=rownames(gadm@data),
                     prod=c(12,26,12,22,0,11,4,5,4,4))
ghaMap <- merge(ghaMap,ghaDF)
ggplot(ghaMap, aes(x=long,y=lat,group=group))+
  geom_polygon(aes(fill=prod))+
  geom_path(color="gray")+
  scale_fill_gradient(low = "light green", high = "dark green")+
  coord_map()