[R-sig-Geo] Colouring maps so that adjacent polygons differ in colour

[Rainer M Krug]

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

On 04/04/11 10:31, Karl Ove Hufthammer wrote:
> 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.

I think the idea is really great. I couldn't try it out yet, but if it
is working, I think it should be added to one of the spatial packages -
maybe sp?

Rainer

> 
> # 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,]
> 


- -- 
Rainer M. Krug, PhD (Conservation Ecology, SUN), MSc (Conservation
Biology, UCT), Dipl. Phys. (Germany)

Centre of Excellence for Invasion Biology
Natural Sciences Building
Office Suite 2039
Stellenbosch University
Main Campus, Merriman Avenue
Stellenbosch
South Africa

Tel:        +33 - (0)9 53 10 27 44
Cell:       +27 - (0)8 39 47 90 42
Fax (SA):   +27 - (0)8 65 16 27 82
Fax (D) :   +49 - (0)3 21 21 25 22 44
Fax (FR):   +33 - (0)9 58 10 27 44
email:      Rainer at krugs.de

Skype:      RMkrug
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.4.10 (GNU/Linux)
Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org/

iEYEARECAAYFAk2bJUIACgkQoYgNqgF2egqM/wCeOe1+UNZaix8PTpLiGY1hZpGD
HHsAn0XGObs/D9QyN8JC3L1V2oFUgMdq
=1EVt
-----END PGP SIGNATURE-----