How to add continuous color legend to an R map made with maps

Question

I'm using the R code shown below, which loads libraries maps and RColorBrewer, to create a map of the world with countries color-coded by population rank. As you can see in the image below, I'm using a green palette in which the darker the green, the larger the population.

I'd like to add a continuous color legend showing the full palette to denote that light green = small population and dark green = large population, but I can't find a way to do it via maps. Could you tell me what is the easiest way to add a continuous color legend (or color key/color scale) to my map?

# Load libraries
library(maps)
library(RColorBrewer)

# Load world data
data(world.cities)

# Calculate world population by country
world.pop = aggregate(x=world.cities$pop, by=list(world.cities$country.etc),
                      FUN=sum)
world.pop = setNames(world.pop, c('Country', 'Population'))

# Create a color palette
palette = colorRampPalette(brewer.pal(n=9, name='Greens'))(nrow(world.pop))

# Sort the colors in the same order as the countries' populations
palette = palette[rank(-world.pop$Population)]

# Draw a map of the world
map(database='world', fill=T, col=palette, bg='light blue')

Answer 1

Within the library SDMTools there is the function legend.gradient

adding this code to the end of your code should give the desired result:

# Draw a map of the world
map(database='world', fill=T, col=palette, bg='light blue')
x = c(-20, -15, -15, -20)
y = c(0, 60, 60, 0)
legend.gradient(cbind(x = x - 150, y = y - 30), 
                  cols = brewer.pal(n=9, name='Greens'), title = "TITLE", limits = "")

You will need to fiddle with the x & y coordinates to get the legend into the desired location however.

EDIT

The x and y coordinates also adjust the shape of the box so I changed the code so that the box shape would not change if you only alter the numbers within the legend.gradient function. Below is what this code should produce

Answer 2

The world map in the maps package is about 30 years old (e.g., has USSR & Yugoslavia).

Plus you have a glitch in your code that causes the overpopulated Greenland that @Jealie noticed (and India is less populated than Antarctica).

You can create a continuousish legend with a modern world using rworldmap.

library(rworldmap)
library(RColorBrewer)

#get a coarse resolution map
sPDF <- getMap()

#using your green colours
mapDevice('x11') #create a map shaped device
numCats <- 100 #set number of categories to use
palette = colorRampPalette(brewer.pal(n=9, name='Greens'))(numCats)
mapCountryData(sPDF, 
               nameColumnToPlot="POP_EST", 
               catMethod="fixedWidth", 
               numCats=numCats, 
               colourPalette=palette)

You can alter the legend adding more labels etc. by doing something like this :

mapParams <- mapCountryData(sPDF, nameColumnToPlot="POP_EST", catMethod="pretty", numCats=100, colourPalette=palette, addLegend=FALSE)

#add a modified legend using the same initial parameters as mapCountryData               
do.call( addMapLegend, c( mapParams
                          , legendLabels="all"
                          , legendWidth=0.5
))

Just briefly to explore the glitch in your code. It occurs because you create a palette for the number of countries in world.cities (239) and then apply it to the number of polygons in the world database from maps (2026). So it probably gets recycled and the colours of your countries have no relation to population. The code below demonstrates the source of your problem.

#find the countries used in the maps world map
mapCountries <- unique( map('world',namesonly=TRUE) )
length(mapCountries)
#[1] 2026
#exclude those containing ':' e.g. "USA:Alaska:Baranof Island"
mapCountries2 <- mapCountries[-grep(':',mapCountries)]
length(mapCountries2)
#[1] 186

#which don't match between the map and world.cities ?
#cityCountries <- unique( world.cities$country.etc )
cityCountries <- world.pop$Country
length(cityCountries)
#[1] 239

#which countries are in the map but not in world.cities ?
mapCountries2[ is.na(match(mapCountries2,cityCountries)) ]
#includes USSR, Yugoslavia & Czechoslovakia