[R-sig-Geo] Colouring maps so that adjacent polygons differ in colour
Karl Ove Hufthammer
karl at huftis.org
Mon Apr 4 10:31:00 CEST 2011
It’s sometimes useful to visualise maps with colours for each
polygon, but where adjacent polygons have different colours.
I recently had to do this, and hope the resulting code may be
of use to others on this list. Note that the code uses a quick
and simple greedy algorithm, which is not guaranteed to find
a colouring using the *fewest* number of colours.
It is *possible* to colour every map using at most four colours
(try searching for ‘four colour theorem’ in your favourite search
engine), but finding the optimal colouring isn’t easy, and an
algorithm implementing this would be much more complicated, and
in some cases rather slow. In practice, the greedy algorithm
usually works very well. Note that sometimes it helps to order
the polygons in some systematic way first, e.g., by the number
of neighbours (the polygon with most neighbours first); see the
example at the end.
I would appreciate any suggestions for improvements to the code.
# Generate colour numbers (indices) for colouring maps
generateMapColours=function(x) {
nb=spdep::poly2nb(x) # Generate neighbours lists
n=length(x) # Number of polygons
cols=numeric(n) # Initial colouring
cols[1]=1 # Let the first polygon have colour 1
cols1n=1:n # Available colour indices
for(i in 2:n)
cols[i]=which.min(cols1n %in% cols[nb[[i]]])
cols
}
## Example
library(rgdal) # Needed for maps in EPSG 4326 not to be interpreted as projected
library(maps) # For map data
library(maptools) # For converting map data to ‘sp’ objects
# Fetch a map
x.map = map("state", fill=TRUE, col="transparent", plot=FALSE)
x.sp = map2SpatialPolygons(x.map, x.map$names, proj4string=CRS("+init=epsg:4326"))
# Plot the map with the calculated colours
cols = generateMapColours(x.sp) # Generate the colour indices
ncols = max(cols) # Number of colours used
library(RColorBrewer) # Nice colour palettes for maps
pal = brewer.pal(ncols, "Set1") # Try ‘display.brewer.all()’ for some alternative palettes
plot(x.sp, col=pal[cols]) # Plot the map
# Ordering the polygons may help.
# For instance, running the colouring algorithm
# on the following produces a four-colouring.
nb = spdep::poly2nb(x.sp)
ord = order(sapply(nb, length), decreasing=TRUE)
x.sp = x.sp[ord,]
--
Karl Ove Hufthammer
More information about the R-sig-Geo
mailing list