## Plot color symbols upon s.class {ade4} outputs

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:

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")) on.exit(par(opar)) 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) box() invisible(match.call()) }

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) { scatterutil.eti.jp<-function (x, y, label, clabel, boxes = TRUE, coul = rep(1, length(x)), horizontal = TRUE) { if (length(label) == 0) return(invisible()) if (is.null(label)) return(invisible()) if (any(label == "")) return(invisible()) 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)) data.frame(x) } opar <- par(mar = par("mar")) par(mar = c(0.1, 0.1, 0.1, 0.1)) on.exit(par(opar)) dfxy <- data.frame(dfxy) if (!is.data.frame(dfxy)) stop("Non convenient selection for dfxy") if (any(is.na(dfxy))) 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) box() invisible(match.call()) }

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

library(ade4) data(deug) 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)

colv<-rep("black",times=length(as.numeric(deug$result))) pchv<-rep(1,times=length(as.numeric(deug$result))) cc<-topo.colors(length(unique(as.numeric(deug$result))))

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

for (j in 1:(length(unique(ccn))) ){ ww<-which(ccn==unique(ccn)[j]) 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)

and superimpose the symbols using s.label.bis

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

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

## Playing with R within a GRASS environment

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. *T. 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.

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.

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 T*omicus piniperda* L. (Coleoptera, Scolytinae). Heredity, 103, 238–247.

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

***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**

mess<-function(X,V,full=TRUE) { messi<-function(p,v){ niinf<-length(which(v<=p)) 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)) return(simi) } E<-extract(x=X,y=1:ncell(X)) r_mess<-X 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<-r_mess[[1]] E<-extract(x=r_mess,y=1:ncell(r_mess[[1]])) 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)`

.

**3-Examples**

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.

library(dismo);library(raster) filename bradypus bradypus files ’/ex’, sep=’’), pattern=’grd’, full.names=TRUE ) predictors predictors<-dropLayer(x=predictors,i=9) reference_points<-extract(predictors,bradypus) mess.out<-mess(X=predictors,V=reference_points,full=TRUE) plot(mess.out)

This is the result:

**4-References**

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. http://CRAN.R-project.org/package=raster

## 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`

.

path<-"/home/rossi/data/blog/data/" map1.n<-paste(path,"map1",".png",sep="") 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`

PlotOnStaticMap(niamey,destfile=map1.n)

Here is is what you get :

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:

newcoord<-LatLon2XY.centered(niamey,lat=13.51204,lon=2.09884,zoom=zoom)

We project the pie chart using these new coordinates:

PlotOnStaticMap(niamey,destfile=map1.n) floating.pie(xpos=newcoord$newX,ypos=newcoord$newY,x=c(20,5), radius=20,col=c("red","blue"))

`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`

.

col<-c("red","blue") for (i in 1:(2)) {col[i]<-col2alpha(color=col[i], alpha = 0.50)} PlotOnStaticMap(niamey,destfile=map1.n) floating.pie(xpos=newcoord$newX,ypos=newcoord$newY,x=c(20,5), radius=20,col=col)

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

Pretty cool he?