[R-sig-Geo] Map colouring via classInt and colour Brewer

[Roger Bivand]

On Wed, 5 Mar 2014, Roger Bivand wrote:

> On Wed, 5 Mar 2014, Alessandra Carioli wrote:
>
>> Dear all,
>> 
>> I have been trying to colour my maps using the same variable at different 
>> points in time employing the classInt and the RColorBrewer palettes. 
>> Although the matching between the map areas and data file is correct (I 
>> have checked), the colouring is wrong every time. I am sure the error is 
>> quite silly but I can’t understand where it lies… I have tried using a 
>> custom made palette or fixed breaks, but the colouring just does not work 
>> the right way! Any help on the matter would be greatly appreciated!
>> 
>> Ale
>> 
>> Library(classInt)
>> #number of class intervals
>> nclassint <- 7
>> # variable to be plotted
>> varofint <- varT1
>
> length(varofint)
>
>
>> #define class intervals
>> cat3 <- classIntervals(varofint, nclassint,style = "jenks")
>> categ <- cat3
>> colpal   = brewer.pal(nclassint,"RdBu")
>
> length(rev(colpal))
>
>
>> # Code for map: I want the reverse of the palette, from blue to red
>> colors <- findColours(categ,rev(colpal))
>> bins <- categ$brks
>> lb <- length(bins)
>> 
>> 
>> plot(shapefile, col=rev(colpal),axes=F)
>
> length(varofint) == length(rev(colpal)) # FALSE, so rev(colpal) will be 
> recycled,
>
> plot(shapefile, col=rep(rev(colpal), length.out=100),axes=F)

This should be:

plot(shapefile, col=rep(rev(colpal), length.out=length(varofint)), ...

I was using North Carolina to replicate with - I wanted to demonstrate the 
consequences of unintended recycling in a reproducible way. The top of the 
code was:

library(maptools)
shapefile <- readShapeSpatial(system.file("shapes/sids.shp",
  package="maptools")[1], IDvar="FIPSNO",
  proj4string=CRS("+proj=longlat +ellps=clrk66"))
varofint <- shapefile$BIR74
library(classInt)
nclassint <- 7
cat3 <- classIntervals(varofint, nclassint, style = "jenks")
categ <- cat3
library(RColorBrewer)
colpal   = brewer.pal(nclassint, "RdBu")
colors <- findColours(categ, rev(colpal))
opar <- par(mfrow=c(2,1))
plot(shapefile, col=rev(colpal), axes=FALSE)
plot(shapefile, col=rep(rev(colpal), length.out=length(varofint)),
  axes=FALSE)
par(opar)

but the intention was:

plot(shapefile, col=colors, axes=FALSE)

Roger

>
> when what you want is:
>
> plot(shapefile, col=colors, axes=F)
>
> Because recycling may be desired, it doesn't generate a warning.
>
> Roger
>
>
>> title("Total Fertility Rate 1981", cex=1)
>> legend(1096133,4844461,fill=rev(colpal),legend=paste(round(bins[-length(bins)],2)
>>                                        ,"-",round(bins[-1],2)),cex=0.6, 
>> bg="white")
>> 
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>> 
>
>

-- 
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no