Archive for the ‘R’ Category

Plot color symbols upon s.class {ade4} outputs

24 August 2012 1 comment

I am used to perform mutlivariate analyses using the R package ade4 by D. Chessel et al. The package features some wonderful graphical functions like s.class or s.label. Unfortunately, these functions not always allow using different colors or pch values when plotting the scatter diagrams with representation of point classes.
The present post deals with that point and will provide some code allowing to get graphical outputs looking like the one below:

A scatter diagram with classes and different symbols & colors

The functions s.label {ade4} does not allow plotting different symbols for the classes and for that reason I have modified the function and transformed it into s.label.bis. It is:

s.label.bis <- function (dfxy, xax = 1, yax = 2, label = row.names(dfxy), clabel = 1,
    pch = 20, cpoint = if (clabel == 0) 1 else 0, boxes = TRUE,
    neig = NULL, cneig = 2, xlim = NULL, ylim = NULL, grid = TRUE,
    addaxes = TRUE, cgrid = 1, include.origin = TRUE, origin = c(0,
        0), sub = "", csub = 1.25, possub = "bottomleft", pixmap = NULL,
    contour = NULL, area = NULL, add.plot = FALSE, col="black", bg=NULL,cex=par("cex"))
    dfxy <- data.frame(dfxy)
    opar <- par(mar = par("mar"))
    par(mar = c(0.1, 0.1, 0.1, 0.1))
    coo <- scatterutil.base(dfxy = dfxy, xax = xax, yax = yax,
        xlim = xlim, ylim = ylim, grid = grid, addaxes = addaxes,
        cgrid = cgrid, include.origin = include.origin, origin = origin,
        sub = sub, csub = csub, possub = possub, pixmap = pixmap,
        contour = contour, area = area, add.plot = add.plot)
    if (!is.null(neig)) {
        if (is.null(class(neig)))
            neig <- NULL
        if (class(neig) != "neig")
            neig <- NULL
        deg <- attr(neig, "degrees")
        if ((length(deg)) != (length(coo$x)))
            neig <- NULL
    if (!is.null(neig)) {
        fun <- function(x, coo) {
            segments(coo$x[x[1]], coo$y[x[1]], coo$x[x[2]], coo$y[x[2]],
                lwd = par("lwd") * cneig)
        apply(unclass(neig), 1, fun, coo = coo)
    if (clabel > 0)
        scatterutil.eti(coo$x, coo$y, label, clabel, boxes)
    if (cpoint > 0 & clabel < 1e-06)
        points(coo$x, coo$y, pch = pch, cex = cex,col=col, bg=bg)

Similarly I changed the ade4 function s.class into s.class.bis which does not plot the “branches” of the stars but only the centroid of each class. It is:

s.class.bis<-function (dfxy, fac, wt = rep(1, length(fac)), xax = 1, yax = 2,
    cstar = 1, cellipse = 1.5, axesell = TRUE, label = levels(fac),
    clabel = 1, cpoint = 1, pch = 20, col = rep(1, length(levels(fac))),
    xlim = NULL, ylim = NULL, grid = TRUE, addaxes = TRUE, origin = c(0,
        0), include.origin = TRUE, sub = "", csub = 1, possub = "bottomleft",
    cgrid = 1, pixmap = NULL, contour = NULL, area = NULL, add.plot = FALSE)
{<-function (x, y, label, clabel, boxes = TRUE, coul = rep(1, length(x)),
    horizontal = TRUE)
    if (length(label) == 0)
    if (is.null(label))
    if (any(label == ""))
    cex0 <- par("cex") * clabel
    for (i in 1:(length(x))) {
        cha <- as.character(label[i])
        cha <- paste(" ", cha, " ", sep = "")
        x1 <- x[i]
        y1 <- y[i]
        xh <- strwidth(cha, cex = cex0)
        yh <- strheight(cha, cex = cex0) * 5/3
        if (!horizontal) {
            tmp <- scatterutil.convrot90(xh, yh)
            xh <- tmp[1]
            yh <- tmp[2]
        if (boxes) {
            rect(x1 - xh/2, y1 - yh/2, x1 + xh/2, y1 + yh/2,
                col = "white", border = coul[i])
        if (horizontal) {
            text(x1, y1, cha, cex = cex0, col = coul[i])
        else {
            text(x1, y1, cha, cex = cex0, col = coul[i], srt = 90)
    f1 <- function(cl) {
        n <- length(cl)
        cl <- as.factor(cl)
        x <- matrix(0, n, length(levels(cl)))
        x[(1:n) + n * (unclass(cl) - 1)] <- 1
        dimnames(x) <- list(names(cl), levels(cl))
    opar <- par(mar = par("mar"))
    par(mar = c(0.1, 0.1, 0.1, 0.1))
    dfxy <- data.frame(dfxy)
    if (!
        stop("Non convenient selection for dfxy")
    if (any(
        stop("NA non implemented")
    if (!is.factor(fac))
        stop("factor expected for fac")
    dfdistri <- f1(fac) * wt
    coul = col
    w1 <- unlist(lapply(dfdistri, sum))
    dfdistri <- t(t(dfdistri)/w1)
    coox <- as.matrix(t(dfdistri)) %*% dfxy[, xax]
    cooy <- as.matrix(t(dfdistri)) %*% dfxy[, yax]
    if (nrow(dfxy) != nrow(dfdistri))
        stop(paste("Non equal row numbers", nrow(dfxy), nrow(dfdistri)))
    coo <- scatterutil.base(dfxy = dfxy, xax = xax, yax = yax,
        xlim = xlim, ylim = ylim, grid = grid, addaxes = addaxes,
        cgrid = cgrid, include.origin = include.origin, origin = origin,
        sub = sub, csub = csub, possub = possub, pixmap = pixmap,
        contour = contour, area = area, add.plot = add.plot)
    if (clabel > 0)
        scatterutil.eti(coox, cooy, label, clabel, coul = col)

Now let look at the example provided with the ade4 function dudi.pca

pca1 <- dudi.pca(deug$tab, scale = FALSE, scan = FALSE)

We now construct different vectors that will be necessary later (this code is necessary because of the data structure of deug)


and use these vectors to create new color and pch vectors that will be used for plotting:

for (j in 1:(length(unique(ccn))) ){
   colv[ww]<-cc[unique(ccn)][j] ; pchv[ww]<-cc2[j]

We now plot the usual scatter plot using s.class

s.class(pca1$li,deug$result,grid=F,cgrid=0,csub=2,cellipse=0,col=cc,cpoint = 0)

output of the s.class function

and superimpose the symbols using s.label.bis

s.label.bis(pca1$li, pch=pchv, clabel = 0, col=colv, add.plot = T)

output of the s.label.bis function

The modified s.class.bis function can be used to redraw the class labels :

s.class.bis(pca1$li, deug$result,csub = 2,cellipse=0, col=cc, cpoint = 0,label=levels(deug$result), add.plot = T)

s.label.bis and s.class.bis can also be used to draw color scatter class graphics such as this:

The MIGCLIM R package

I just saw the paper by Engler et al (2012) published a few days ago presenting the R version of MIGCLIM. The MIGCLIM R package allows the integration of dispersal constraints into projections of species distribution models. This is a very promising tool and I will test it as soon as possible.

Robin Engler, Wim Hordijk, Antoine Guisan (2012) The MIGCLIM R package – seamless integration of dispersal constraints into projections of species distribution models. Ecography. In press. DOI: 10.1111/j.1600-0587.2012.07608.x Article first published online: 3 AUG 2012

paper abstract

Categories: R Tags: ,

Playing with R within a GRASS environment

14 August 2012 1 comment

Although there are some powerful GIS utilities in R, I prefer using GRASS to manage my GIS data while I use R to perform scientific computing.

In fact GRASS and R can be very simply interfaced by means of the R package spgrass6.

Under linux operating system (actually Kubuntu Natty Narwhal) I open the console, launch GRASS by typing grass and then select a location (see GRASS manual for details).


Once GRASS is running, R can be launched from within the GRASS console…

From now on I am working in R and I can use the library spgrass6 in order e.g. to read or write rasters into the GRASS system.

As an example we will consider the case of the pine shoot beetle Tomicus piniperda (Coleoptera: Curculionidae: Scolytinae). Readers can find the whole story in Horn et al. (2012) available here. Tomicus piniperdaT. piniperda is present throughout Europe. Interestingly, it has long been assumed to be present in North Africa too although molecular studies have recently shown that T. piniperda only rarely occurs in these regions (Horn et al., 2006, 2009).

We have a set of localities where T. piniperda was recorded as present or, on the contrary, has been searched and was absent (true absence). GRASS is used to host the data and produce the graphical outputs.

Sites where T. piniperda was either present (black circle) or absent (open circle).

Starting R from the GRASS console and using the dismo package allows us to fit a SDM (Species Distribution Model) and build a raster with the probabilities of presence e.g. using the function predict. This raster layer can be expressed as presence/absence and the resulting data written within the GRASS system using the function writeRAST6 from package spgrass6.
The raster can now be plotted using the GRASS interface and its utilities.

GRASS window showing the occurences of T. piniperda and the associated predicted distribution derived from a SDM run in R

Horn, A., Kerdelhué, C., Lieutier, F., Rossi, J.-P. (2012). Predicting the distribution of the two bark beetles Tomicus destruens and Tomicus piniperda in Europe and the Mediterranean region. Agricultural and Forest Entomology, in press, DOI: 10.1111/j.1461-9563.2012.00576.x.

Horn, A., Roux-Morabito, G., Lieutier, F. & Kerdelhué, C. (2006) Phylogeographic structure and past history of the circum-Mediterranean species Tomicus destruens Woll. (Coleoptera: Scolytinae). Molecular Ecology, 15, 1603–1615.

Horn, A., Stauffer, C., Lieutier, F. & Kerdelhué, C. (2009) Complex postglacial history of the temperate bark beetle Tomicus piniperda L. (Coleoptera, Scolytinae). Heredity, 103, 238–247.

Computing the “Multivariate Environmental Similarity Surfaces” (MESS) index in R

13 August 2012 9 comments

*edit oct 03 2012 *
The mess function is now part of the dismo package developped by Robert J. Hijmans, Steven Phillips, John Leathwick and Jane Elith (dismo 0.7-23)***

This post is about  the computation of the “Multivariate Environmental Similarity Surfaces” aka MESS in R. The MESS index was proposed by Elith et al (2010) [Methods in Ecology & Evolution 2010, 1, 330–342]. MESS can be computed with the Maxent software but I  wanted to perform all my data treatments within R/GRASS environment. For that reason I wrote a R function called mess.R.

Note that the function requires the R package raster developed by R. Hijmans and J. van Etten. Both the function and a tutorial can be downloaded for my site.

The present post provides some explanations, the mess.R function and some examples.

1- What is MESS?
MESS stands for Multivariate Environmental Similarity Surfaces. It is an index of similarity reporting the closeness of a point described by a set of environmental attributes to the distribution of these attributes within a population of reference points. The MESS approach has been proposed by Elith et al (2010). It works as BIOCLIM does but provides negative values which allows one to differentiate levels of dissimilarity. This is particularly valuable when considering points associated to values lying outside the range of the reference points.

2-The mess.R function

    f    if(f==0) simi<-100*(p-min(v))/(max(v)-min(v))
    if(0    if(50<=f & f<100) simi<-2*(100-f)
    if(f==100) simi<-100*(max(v)-p)/(max(v)-min(v))
for (i in 1:(dim(E)[2])) {
    e<-data.frame(E[,i]) ; v<-V[,i]
    r_mess[[i]][]<-apply(X=e, MARGIN=1, FUN=messi, v=v)
rmess[]<-apply(X=E, MARGIN=1, FUN=min)
if(full==TRUE) {
    out     layerNames(out)<-c(layerNames(X),"mess")
if(full==FALSE) out return(out)

The function mess.R allows the computation of the MESS index for a set of raster objets. raster objects are defined within the package raster developped by Hijmans & van Etten (2011). mess.R therefore requires that the package raster is properly installed and loaded using the command library(raster).

The following example uses the Bradypus data hosted in the package dismo (Hijmans et al., 2012). As with raster, dismo must be installed and loaded.

filename bradypus bradypus files ’/ex’, sep=’’), pattern=’grd’, full.names=TRUE )
predictors predictors<-dropLayer(x=predictors,i=9)

This is the result:

Elith J., Kearney M., & Phillips S. 2010. The art of modelling range-shifting species. Met. Ecol. Evol 1 :330-342.

Hijmans R.J. & van Etten J. (2011). raster : Geographic analysis and modeling with raster data. R package version 1.9-58.

Hijmans R.J., Phillips S., Leathwick J. & Elith J. (2012). dismo : Species distribution modeling. R package version 0.7-17.

Categories: R Tags: ,

Display pie charts with varying alpha transparencies upon a Google static map

It is sometimes interesting to display pie charts on a map showing geographic details such as Google static maps. This post will illustrate :

  • how to download and manipulate Google static maps using the R package {RgoogleMaps} by M. Loecher
  • how to overlay floating pie charts upon a map using the R package {plotrix} by J. Lemon
  • how to change the alpha value associated to colors {seqinr} by D. Charif, and J.R. Lobry

As an example we will focus on the town of Niamey where my colleagues Madougou Garba and Gauthier Dobigny (Centre Régional AGRHYMET, Département Formation Recherche, Niamey) are currently studying urban rodents communities. Lets say we want to display pie charts corresponding to sampling locations where the rodent community has been sampled using traps.

First we download the static map from Google server using the function GetMap from the package RgoogleMaps.

zoom<-12 ; lonc<-2.114668 ; latc<-13.520508 ; center<-c(latc,lonc) ; maptype<-"satellite"
niamey <- GetMap(center=center, zoom=zoom, maptype=maptype,destfile=map1.n)

Note that lonc and latc are the coordinates of a points roughtly located at the center of the town. It indicates where the center of the static map should be located.

Once the data are available through the R object niamey we can display it on sreeen using PlotOnStaticMap


Here is is what you get :

Google static map of Niamey

Now we come to the pie charts. Lets say we have a sample of the urban rodent community collected in a location which coordinates are latitude=13.51204 and longitude=2.09884. Imagine that 20 individuals of R. rattus and 5 individuals of M. natalensis were recorded at that location.

We can plot a point that shows the spatial location of the locality

PlotOnStaticMap(niamey,lat=13.51204, lon=2.09884, pch=3,cex = 3,col="red",verbose=0,destfile=map1.n)


If we want a pie chart showing the proportion of each rodent species we use the function floating.pie(xpos,ypos,x,...) where xpos and ypos are the coordintates of the centre of the pie and x is a vector with rodent species counts.
The trick is that xpos and ypos should correspond to the map tile coordinates system and cannot be given in raw degrees.
The package RgoogleMaps provides the function LatLon2XY.centered allowing the conversion:


We project the pie chart using these new coordinates:


radius controls the size of the pie chart, and it can be adjusted to plot a set of pie chart with size proportional to the sample size for example.

When a large number of pie charts are plotted on the same map they may be superimposed. In that case having more or less transparent colors (i.e. changing the alpha) may be interesting.
This is achieved using the function col2alpha from the package seqinr.

for (i in 1:(2)) {col[i]<-col2alpha(color=col[i], alpha = 0.50)}

Here is the final result, note that the river Niger can now be seen through the pie chart:

Pretty cool he?

Categories: R Tags: ,