Clarke and Ainsworth's BIOENV and BVSTEP (and BIO-BIO etc...)

Nonmetric Multidimensional Scaling (NMDS) plot of vegetation sample dissimilarities with best correlating environmental variables (left) and species (right) plotted as vectors (datasets "varespec" and "varechem" from the package vegan)

The R package "vegan" contains a version of Clarke and Ainsworth's (1993) BIOENV analysis allowing for the comparison of distance/similarity matrices between two sets of data having either samples or variables in common. The typical setup is in the exploration of environmental variables that best correlate to sample similarities of the biological community (e.g. species biomass or abundance). In this case, the similarity matrix of the community is fixed, while subsets of the environmental variables are used in the calculation of the environmental similarity matrix. A correlation coefficient (typically Spearman rank correlation coefficient) is then calculated between the two matrices and the best subset of environmental variables can then be identified and further subjected to a permutation test to determine significance.

This can be a very helpful analysis in the exploration of the often highly dimensional space of community samples. The method is also widely accepted by the scientific community due to its flexibility across a wide variety of data and is completely non-parametric - Clarke and Ainsworth's (1993) paper describing the method has 674 citations on Google Scholar at the time of this posting. 

The R package "vegan" incorporates this routine in the function bioenv(). An example of a BIOENV exploration between vegetation community data (dataset "varespec" in the vegan package) and the environmental data (dataset "varechem" in the vegan package) :

Whales, plankton migrate across Northwest Passage

Phytoplankton species not seen in the Atlantic for at least 800,000 years are turning up in recent years due to an increase in water transport through the Northwest Passage. While this is not yet cause for alarm, increases in species invasions from one oceanic basin to another could cause larger shifts in marine ecosystems over time - read more here

Image color palette replacement

Here is an example of a function I wrote to change the color palette used in an image. The above example comes from a black and white original, although color images can also be used. The function first converts the image to grayscale in order to have levels of color intensity between 0-255. Using a new color palette with 256 color levels, the gray levels are replaced with a rgb (red, blue, green) vector from the new palette. The results can be very strange...
The package biOps is required for reading and writing the .jpeg files.

the function...

Simulating CMYK mis-registration printing

I recently came across a poster advertising a children's production of Shakespeare's The Tempest where they purposely used an effect to mimic a mis-registration in CMYK printing. You have probably seen this before as a slight offset in one of the 4 colors (cyan, magenta, yellow, and black).
     The CMYK color model is "subtractive" in that you end up with a white color (when the paper color is white) when no colors are printed. The opposite is an "additive" color model, such as RGB (red, green, blue), which results in black when none of these three color channels are added. This is more typically used in imaging on lit screens (e.g. color creation in R using the rgb() function).
     I wanted to try simulating this type of mis-registration in R and came up with the following function. For images with white backgrounds, the "subtractive" shift will look best while an "additive" shift works best for black backgrounds. The results are essentially the same, but you can eliminate some color channel striping on the image borders by choosing one or the other.
     This is probably much easier to do in a photo editting program, but I had fun with it nonetheless. I used the excellent package biOps for some of its image reading and manipulation functions.

...the function

map.xyz(): interpolation of XYZ data and projection onto a map

     I am still struggling to get a grasp of R's mapping capabilities. Part of my frustration lies in the fact that I often work on areas near the poles, which complicates interpolation across the 180 degree line. For smaller areas, interpolation can be done using the interp() function in the package akima. I have taken the results from interp and projected the image onto a map. You will need the akima, maps, and mapproj packages and the functions new.lon.lat(), earth.dist(), and pos2coord().


As an example I have mapped the distance from Mecca, Saudi Arabia:

The function...

Array position to matrix coordinates conversion

#A function that is sometimes useful in determining the 
#coordinate(i.e. row and column number) of a matrix position
#(and vice-versa). 
#Either a vector of positions ("pos") 
#OR a 2 column matrix of matrix coordinates, ("coord", i.e. cbind(row,col)), 
#AND the matrix dimentions must be supplied (dim.mat, i.e. c(nrow,ncol)).
pos2coord<-function(pos=NULL, coord=NULL, dim.mat=NULL){
 if(is.null(pos) & is.null(coord) | is.null(dim.mat)){
  stop("must supply either 'pos' or 'coord', and 'dim.mat'")
 }
 if(is.null(pos) & !is.null(coord) & !is.null(dim.mat)){
  pos <- ((coord[,2]-1)*dim.mat[1])+coord[,1] 
  return(pos)
 }
 if(!is.null(pos) & is.null(coord) & !is.null(dim.mat)){
  coord <- matrix(NA, nrow=length(pos), ncol=2)
  coord[,1] <- ((pos-1) %% dim.mat[1]) +1
  coord[,2] <- ((pos-1) %/% dim.mat[1]) +1
  return(coord)
 }
}

Created by Pretty R at inside-R.org

R functions for Earth geographic coordinate calculations

Here are some functions that I regularly use for geographic data (e.g. binning, filtering, calculation of new positions etc.).

#distance in kilometers between two long/lat positions (from "fossil" package)
earth.dist <- function (long1, lat1, long2, lat2) 
{
    rad <- pi/180
    a1 <- lat1 * rad
    a2 <- long1 * rad
    b1 <- lat2 * rad
    b2 <- long2 * rad
    dlon <- b2 - a2
    dlat <- b1 - a1
    a <- (sin(dlat/2))^2 + cos(a1) * cos(b1) * (sin(dlon/2))^2
    c <- 2 * atan2(sqrt(a), sqrt(1 - a))
    R <- 6378.145
    d <- R * c
    return(d)
}

Created by Pretty R at inside-R.org

#degree bearing between two long/lat positions (from "fossil" package)
earth.bear <- function (long1, lat1, long2, lat2) 
{
    rad <- pi/180
    a1 <- lat1 * rad
    a2 <- long1 * rad
    b1 <- lat2 * rad
    b2 <- long2 * rad
    dlon <- b2 - a2
    bear <- atan2(sin(dlon) * cos(b1), cos(a1) * sin(b1) - sin(a1) * 
        cos(b1) * cos(dlon))
    deg <- (bear%%(2 * pi)) * (180/pi)
    return(deg)
}

Created by Pretty R at inside-R.org

new.lon.lat <-
function (lon, lat, bearing, distance) 
{
    rad <- pi/180
    a1 <- lat * rad
    a2 <- lon * rad
    tc <- bearing * rad
    d <- distance/6378.145
    nlat <- asin(sin(a1) * cos(d) + cos(a1) * sin(d) * cos(tc))
    dlon <- atan2(sin(tc) * sin(d) * cos(a1), cos(d) - sin(a1) * 
        sin(nlat))
    nlon <- ((a2 + dlon + pi)%%(2 * pi)) - pi
    npts <- cbind(nlon/rad, nlat/rad)
    return(npts)
}

Created by Pretty R at inside-R.org

#tells which lon lat positions are within the defined limits to the west, east, north, and south
lon.lat.filter <-
function (lon_vector, lat_vector, west, east, north, south) 
{
 if(west>east) {
  lon_vector_new=replace(lon_vector, which(lon_vector<0), lon_vector[which(lon_vector<0)]+360)
  east_new=east+360
 } else {
  lon_vector_new=lon_vector
  east_new=east
 }
  hits=which(lon_vector_new < east_new & lon_vector_new > west & lat_vector < north & lat_vector > south)
 return(hits)
} 

Created by Pretty R at inside-R.org