GIS, an acronym that brings joy to some and strikes fear in the heart of those not interested in buying expensive software. Luckily fight or flight can be saved for another day because you don’t need to be a GIS jock with a wad of cash to work with spatial data and make beautiful plots. A computer and internet connection should be all you need. In this post I will show how to

• Get a machine ready to use R to work with geospatial data
• Describe what type of data can be used and some of the exciting sources of free GIS data
• Use data from the Washington Department of Natural Resources to generate some pretty plots by working through an example script

Getting your Machine Ready

First, if you do not have R on your computer you can download it here. Even better, download Rstudio, an incredibly efficient and easy to use interface for working with R available here.  Working with R studio is highly recommended and will be more clearly outlined in this post. R does not support working with spatial data straight out of the box so there are a couple of packages that need to be downloaded to get R working with spatial data. The two packages required are ‘sp’ and ‘rgdal’.  We will also use a third package, ‘rgeos’ for some fancy geospatial tricks. Unfortunately, the latest release of the sp package is not compatible with the latest version of R — v 3.0 at this time. When the Rstudio install prompts to install R, download version 2.15.3.

First intall and open RStudio. To add the required packages open RStudio and click on the “packages” tab in the lower right hand panel of the interface. The lower right window will show a list of packages that come with a standard download of RStudio. Some of the packages will have check marks next to them, this means that those libraries are loaded and ready to be used. If you just downloaded R for the first time sp and rdgal will not be on that list, click on the “Install Packages” button. Make sure the “Install from” option is set to “Repository (CRAN)” and type “sp” into  the “Packages” space. Check the “Install Dependencies” option and download! By checking the “Install Dependencies” option packages sp needs to function properly will automatically be downloaded. Download rgdal in the same way and you have the tools needed to start!  Download rgeos as well and you can run the portion of the example script that uses centroids.

Sp and Rgdal abilities and sources of Data

Rgdal is what allows R to understand the structure of shapefiles by providing functions to read and convert spatial data into easy-to-work-with R dataframes. Sp enables transformations and projections of the data and provides functions for working with the loaded spatial polygon objects. The take away should be that these are powerful libraries that allow R to use .shp files. The US Geological Survey, the National Park Serivce, and the Washington State Department of Natural Resources are a just a few examples of organizations that make enormous stockpiles of spatial data available to the public.

Example Script

The following code uses Watershed Resource Inventory Area (WRIA) spatial data from the Washington State Department of Ecology.  This dataset contains information on Washington State’s different water resource management areas.  Exactly what information is stored in the shapefiles will be explored using R! If a function or any of the code looks mysterious try  ”? mysteriousFunctionName()” and handy documentation will fill you in on what a function does. Lets start using R to investigate the data.  Just cut and paste the following code into RStudio.

# WRIA example by Steven Brey @ mazamascience.com

# load the required libraries 
library(sp) 
library(rgdal)

# First load the data from the Washington department of ecology website

# Data source
#
# Data:      ftp://www.ecy.wa.gov/gis_a/hydro/wria.zip
# Metadata:  http://www.ecy.wa.gov/services/gis/data/hydro/wria.htm

# Create a directory for the data
localDir <- 'R_GIS_data'
if (!file.exists(localDir)) {
  dir.create(localDir)
}
# Download the 5mb file of WRIA data
url <- 'ftp://www.ecy.wa.gov/gis_a/hydro/wria.zip'
file <- paste(localDir,basename(url),sep='/')
if (!file.exists(file)) {
  download.file(url, file)
  unzip(file,exdir=localDir)
}
# Show the unzipped files 
list.files(localDir)

# layerName is the name of the unzipped shapefile without file type extensions 
layerName <- "WRIA_poly"  
# Read in the data
data_projected <- readOGR(dsn=localDir, layer=layerName) 

# What is this thing and what's in it?
class(data_projected)
slotNames(data_projected)
# It's an S4 "SpatialPolygonsDataFrame" object with the following slots:
# [1] "data"        "polygons"    "plotOrder"   "bbox"        "proj4string"

# What does the data look like with the default plotting command? 
plot(data_projected)

We don’t know yet what the different ‘slots’ contain but Holy Smokes!  That plot looks like Washington. The plot is not pretty (yet) but in just a few lines of code you can already make a primitive map of the polygons in this dataset. The fact that this object is a ~DataFrame suggests that we can treat it like an R dataframe.  Lets find out what data it has, select the interesting columns and work with only those.  (Not surprisingly, the @data slot is where the actual dataframe exists but many dataframe methods also work on the entire object.)

We’ll go to the trouble of renaming variables, reprojecting the data to “longlat” and then saving the data to a .RData file for easy loading in the future.

# Could use names(data_projected@data) or just:
names(data_projected)
#["WRIA_ID"    "WRIA_NR"    "WRIA_AREA_" "WRIA_NM"    "Shape_Leng" "Shape_Area"

# Identify the attributes to keep and associate new names with them
attributes <- c("WRIA_NR", "WRIA_AREA_", "WRIA_NM")
# user friendly names 
newNames <- c( "number", "area", "name")

# Subset the full dataset extracting only the desired attributes
data_projected_subset <- data_projected[,attributes]

# Assign the new attribute names
names(data_projected_subset) <- newNames

# Create a dataframe name (potentially different from layerName)
data_name <- "WRIA"

# Reproject the data onto a "longlat" projection and assign it to the new name
assign(data_name,spTransform(data_projected_subset, CRS("+proj=longlat")))

# NOTE: If using assign() above gave you an error it is likely the version of 
# NOTE: R you are using does not currently support the sp package. Use R 
# NOTE: 2.15.3 instead. 

# The WRIA dataset is now projected in latitude longitude coordinates as a
# SpatialPolygonsDataFrame.  We save the converted data as .RData for faster
# loading in the future.
save(list=c(data_name),file=paste(localDir,"WAWRIAs.RData",sep="/"))

# Upon inspecting the metadata you can see that the first 19 areas in WRIA
# surround Puget Sound. The names of the first 19 watersheds in WRIA are
WRIA$name[1:19]

# For fun, save a subset including only only these 19 areas
PSWRIANumbers <- c(1:19)
WRIAPugetSound <- WRIA[WRIA$number %in% PSWRIANumbers,]
# Sanity check, plot WRIAPugetSound to make sure it looks like the subset we want
plot(WRIAPugetSound)

# Save Puget Sound data as an RData file which we can quickly "load()"
data_name <- "WRIAPugetSound"
save(list=c(data_name),file=paste(localDir,"WRIAPugetSound.RData",sep="/"))

We have now saved both WRIA and WRIAPugetSound data as R spatial objects that can easily be loaded! Now the real fun begins, those of you who have been waiting for pretty plots are about to be rewarded. The rest of the script is a walk through of some of the fun analysis and basic figure creating that can easily be done with the converted WRIA data.

# Load the data that was converted to SpatialPolygonsDataFrame
# NOTE: This can be skipped (but does not have to be) if the spatial
# NOTE: objects are still in your workspace.

# Load and show the names of the attributes in WRIA
file <- paste(localDir,"WAWRIAs.RData",sep="/")
load(file) 
names(WRIA)

file <- paste(localDir,"WAWRIAs.RData",sep="/")
load(file)
names(WRIAPugetSound)

# Sweet.  We can see that this is the WRIA dataset we saved earlier

# NOTE: For more advanced users, slotNames(WRIA) will list the structures 
# in WRIA. Using the @ command allows you to grab a particular slot from the
# spatial object.  If you really want the HUGE AND GORY details of what's
# in this object, you can examine the full structure with str(WRIA).

# Here is how you would extract the contents of slots for easier use.
WriaData <- WRIA@data
WriaBbox <- WRIA@bbox

# We have a pretty good idea of what kind of data we are working with 
# and what it looks like. Now its time for the data to answer some
# questions and tell us a story.

# What is the biggest water resource area in Washington? 
maxArea <- max(WRIA$area)
# Create a 'mask' identifying the biggest area so we can find out its name
# NOTE:  Eaach 'mask' we create is a vector of TRUE/FALSE values that we
# NOTE:  will use to subset the dataframe.
biggestAreaMask <- which(WRIA$area == maxArea)
biggestAreaName <- WRIA$name[biggestAreaMask]
biggestAreaName

# Create a SpatialPolygonsDataFrame subset
biggestArea <- WRIA[biggestAreaMask,]

# plot biggest area in blue with a title
plot(biggestArea, col="blue")
title(biggestAreaName) 

# NOTE: Many more plot arguments can be explored by investigating 
# NOTE: the "SpatialPolygons" "plot-method" in the sp package

# I have heard of a water resource management area in Washington State
# called Pend Oreille.  Where is it located in this dataframe?
which(WriaData$name == "Pend Oreille")

# Now we have isolated the watershed with the largest area as well as the
# fabled Pend Oreille.  Lets figure out how to highlight these regions when
# plotting all  regions. I have also heard that Lake Chelan is Beautiful.
# Lets isolate it as well.

# Each of the following makes a spatialPolygonsDataFrame subset, selecting 
# a specific region based on some selected attribute in WRIA.

WRIA_puget <- WRIA[WRIA$number %in% PSWRIANumbers,]
WRIA_chelan <- WRIA[WRIA$name == "Chelan",]
WRIA_Pend_Oreille <- WRIA[WRIA$name == "Pend Oreille",]

# Check out what they look like plotted individually 
plot(WRIA_puget)
plot(WRIA_chelan)
plot(WRIA_Pend_Oreille)

# For fun we will make 8 different watersheds 8 different colors!
watersheds <- c(1:8)
watershed_colors <- c("burlywood1","forestgreen","burlywood3","darkolivegreen3",
                      "cadetblue4","sienna3","cadetblue3","darkkhaki")

watershed_color_variation <- WRIA[WRIA$number %in% watersheds,]

# Plot some of the created spatial objects together
plot(WRIA)
plot(WRIA_puget,add=TRUE,col="navy")
plot(WRIA_chelan,add=TRUE,col="red2")
plot(watershed_color_variation, add=TRUE, col=watershed_colors)
plot(WRIA_Pend_Oreille,add=TRUE,col="red4")

# NOTE:  gCentroid is from the 'rgeos' package
library(rgeos)
# Find the center of each region and label lat and lon of centers
centroids <- gCentroid(WRIA, byid=TRUE)
centroidLons <- coordinates(centroids)[,1]
centroidLats <- coordinates(centroids)[,2]

# Find regions with center East of -121 longitude (East of Cascade mountains) 
eastSideMask <- centroidLons > -121

# Create spatialPolygonsDataFrame for these regions
WRIA_nonPacific <- WRIA[eastSideMask,]

# Find watersheds with area 75th percentile or greater
basinAreaThirdQuartile <- quantile(WRIA$area,0.75)
largestBasinsMask <- WRIA$area >= basinAreaThirdQuartile

WRIA_largest <- WRIA[largestBasinsMask,]

# To get legend and labels to fit on the figure we can change the size of the
# area plotting. bbox(WRIA) shows the bounding lat and long of WRIA.
bBox <- bbox(WRIA)
ynudge <- 0.5
xnudge <- 3
xlim <- c(bBox[1,1] + xnudge , bBox[1,2] - xnudge)
ylim <- c(bBox[2,1] - ynudge, bBox[2,2] )

# Plot some of the different spatial objects and show lat-long axis, 
# label watersheds
plot(WRIA,axes=TRUE,xlim=xlim,ylim=ylim)
plot(WRIA_nonPacific,add=TRUE,col="red") 
plot(WRIA_puget,add=TRUE,col="navy")
plot(WRIA_largest,add=TRUE,density=20,col="green1")
text(centroidLons, centroidLats, labels=WRIA$number, col="blue", cex=.7)

title(main="Washington Watersheds")
labels <- c("Puget Sound Watersheds", "Washington's biggest watersheds",
            "drain to Pacific via Columbia river") 
label_cols <- c("navy","green1","red")
legend("bottomleft", NULL, labels, fill=label_cols, density, bty="n", cex=.8)

Who needs expensive GIS software?

Another useful package for plotting spatial data is the “maptools” library available from from the Cran repository. The code below requires the maptools package so make sure it is installed before running the code.

# Lets get a taste of what some of the built in plotting functions in the map
# tools package can do

# First load
library("maptools")

allWashingtonAreas <- spplot(WRIA, zcol="area", identify=TRUE)

# by saving as object more than one plot can be shown on a panel
pugetSoundAreas <- spplot(WRIA_puget, zcol="area", identify=TRUE)

# Open new plot frame
frame()
print(allWashingtonAreas, position=c(0,.5,.5,1), more=T)
print(pugetSoundAreas, position=c(.5,.5,1,1), more=T)

These examples only scratch the surface of the plotting and analysis potential available in R using the sp, rgdal, and maptools packages.  Exploring the plot-methods in sp is highly recommended before investing a large amount of time into scripting a plot on your own.

Congratulations! You now know how to find out what is in geospatial data and how to make plots using R.

Happy Mapping!

1. Write programs for people, not computers.
   1. a program should not require its readers to hold more than a handful of facts in memory at once
   2. names should be consistent, distinctive and meaningful
   3. code style and formatting should be consistent
   4. all aspects of software development should be broken down into tasks roughly an hour long
2. Automate repetitive tasks.
   1. rely on the computer to repeat tasks
   2. save recent commands in a file for re-use
   3. use a build tool to automate scientific workflows
3. Use the computer to record history.
   1. software tools should be  used to track computational work automatically
4. Make incremental changes.
   1. work in small steps with frequent feedback and course correction
5. Use version control.
   1. use a version control system
   2. everything that has been created manually should be put in version control
6. Don’t repeat yourself (or others).
   1. every piece of data must have a single authoritative representation in the system
   2. code should be modularized rather than copied and pasted
   3. re-use code instead of rewriting it
7. Plan for mistakes.
   1. add assertions to programs to check their operation
   2. use an off-the-shelf unit testing library
   3. use all available oracles when testing programs
   4. turn bugs into test cases
   5. use a symbolic debugger
8. Optimize software only after it works correctly.
   1. use a profiler to identify bottlenecks
   2. write code in the highest-level language possible
9. Document design and purpose, not mechanics.
   1. document interfaces and reasons, not implementations
   2. refactor code instead of explaining how it works
   3. embed the documentation for a piece of software in that software
10. Collaborate.
    1. use pre-merge code reviews
    2. use pair programming when bringing someone new up to speed and when tackling particularly tricky problems

The only extra I would have included would be:

11. Maintain and update older code.

If you are still hesitant to go to the original article, go there for the 67 references to other books and articles that discuss scientific computing.  Like I said, this article is a “Best Ever”.

The post titled Installing Packages described the basics of package installation with R.  The process is wonderfully simple when everything goes well.  But it can be maddening when it does not.  Error messages give a hint as to what went wrong but do not necessarily tell you how to resolve the problem.  This post will collect some of the error messages we’ve encountered while installing R packages and describe the reasons for the error and the workarounds we’ve found.

1) Older version of R

Warning message:
In install.packages(c("sp")) : package ‘sp’ is not available

This is the message that you get when the CRAN package you’re interested in requires a more recent version of R than you have.  Remember, the default behavior of install.packages() is to grab the latest version of a package.

In this case you have to poke around in the “Old sources” link on the CRAN page for that package and use trial-and-error to find an older version of the package that will work with your version of R.  You should start by determining what version of R you have:

$ R --version
R version 2.8.1 (2008-12-22)

This version of R was released at the end of 2008 and any version of the “sp” package released in 2008 should work. At least some of the 2009 releases should also work. Perusing the sp archive, we might try installing version 0.9-37, the last of the 0.9-3x series which was released in May of 2009:

$ wget http://cran.r-project.org/src/contrib/Archive/sp/sp_0.9-37.tar.gz
$ sudo CMD INSTALL sp_0.9-37.tar.gz
...
$ # Success!

2) Unable to execute files in /tmp directory

ERROR: 'configure' exists but is not executable -- see the 'R Installation and Administration Manual'

By default, R uses the /tmp directory to install packages.  On security conscious machines, the /tmp directory is often marked as “noexec” in the /etc/fstab file.  This means that no file under /tmp can ever be executed.  Packages that require compilation or that have self-inflating data will fail with the error above.  One such package is RJSONIO.

The solution is to set the TMPDIR environment variable which R will use as the compilation directory. For csh shell:

$ mkdir ~/tmp
$ setenv TMPDIR ~/tmp

And for bash:

$ mkdir ~/tmp
$ export TMPDIR=~/tmp

Other problems

Please leave comments describing other package installation problems and solutions you’ve encountered to help us build a more complete listing.

Yes, this post condenses 50+ hours of learning into a 15 minute tutorial.  Read ‘em and weep.  (That is, you read while I weep.)

OK.  For the last week I’ve been learning how to call C code as documented in previous posts:  Calling C code “Hello World!”, .Call(“hello”) and Calling C code with Rcpp.  Now comes that task that puts this hard won knowledge to use — providing an R interface to a library of C code as part of a package.  You are encouraged to have Writing R Extensions open in another window and use search whenever you see something you want more information about.  In an effort to keep this exercise to 15 minutes we will stick to the task at hand and not do a whole lot of explaining.  It is assumed that you already have the required compilers installed.

1) Download and test your C library

We’ll be working with the “libmseed” library that reads binary seismic data.  Download the .gz file from the IRIS DMC Softare Library and then compile and test it with the example code.  (Read the README files on your own time.  Right now we have work to do.)

$ cd ~/Downloads
$ wget http://www.iris.edu/pub/programs/libmseed-2.7.tar.gz
--2012-11-16 16:18:11--  http://www.iris.edu/pub/programs/libmseed-2.7.tar.gz
Resolving www.iris.edu... 128.95.166.129
Connecting to www.iris.edu|128.95.166.129|:80... connected.
HTTP request sent, awaiting response... 200 OK
Length: 166186 (162K) [application/x-gzip]
Saving to: “libmseed-2.7.tar.gz”
...
$ tar -xzf libmseed-2.7.tar.gz 
$ cd libmseed
$ make
cc -c -o fileutils.o fileutils.c
cc -c -o genutils.o genutils.c
cc -c -o gswap.o gswap.c
cc -c -o lmplatform.o lmplatform.c
cc -c -o lookup.o lookup.c
cc -c -o msrutils.o msrutils.c
cc -c -o pack.o pack.c
cc -c -o packdata.o packdata.c
cc -c -o traceutils.o traceutils.c
cc -c -o tracelist.o tracelist.c
cc -c -o parseutils.o parseutils.c
cc -c -o unpack.o unpack.c
cc -c -o unpackdata.o unpackdata.c
cc -c -o selection.o selection.c
cc -c -o logging.o logging.c
rm -f libmseed.a
ar -csq libmseed.a fileutils.o genutils.o gswap.o lmplatform.o lookup.o msrutils.o pack.o packdata.o 
traceutils.o tracelist.o parseutils.o unpack.o unpackdata.o selection.o logging.o
$ cd example
$ make
cc -I.. -c msview.c
cc -I.. -o msview msview.o -L.. -lmseed
cc -I.. -c msrepack.c
cc -I.. -o msrepack msrepack.o -L.. -lmseed
$ ./msview test.mseed
IU_COLA_00_LHZ, 000001, M, 512, 112 samples, 1 Hz, 2010,058,06:50:00.069539
...
IU_COLA_00_LHZ, 000001, M, 512, 27 samples, 1 Hz, 2010,058,07:59:33.069538
$

[02:00 minutes]

Congratulations are in order.  (No, not for you!  For the smart folks who wrote libmseed!) ;-)

2) Set up your package structure

Yes, we’re going to do that now.  You could take incremental steps and practice writing an R callable version of example/msview.c and then try to compile it with R CMD SHLIB, then give up and write a Makefile that contains all of the R CMD SHLIB flags.  But I’m here to tell you that most of that will happen auto-magically if you set up a package first and let R worry about the details.

Here is the minimal package structure you need.  We’ll create a new directory named “mseed” and copy all of the “libmseed” code into it:

$ cd ~
$ mkdir mseed
$ cd mseed
$ mkdir demo
$ mkdir inst
$ mkdir inst/extdata
$ mkdir R
$ mkdir src
$ cp -r ~/Downloads/libmseed ./src

[02:40]

(Curiosity might lead a person to search through your already open copy of “Writing R Extensions” to figure out what “extdata” is used for.  In the interest of time, please do that later.)

3) Add Makevars to your mseed/src/ directory

You need to educate R that you are going to compile code that depends upon the libmseed.a library.  The easiest way to do this is by creating a file named mseed/src/Makevars::

# See Section 1.2.1 "Using 'Makevars'" of Writing R Extensions
# cran.r-project.org/doc/manuals/R-exts.pdf

PKG_CFLAGS=
PKG_CPPFLAGS=-Ilibmseed
PKG_LIBS=-Llibmseed -lmseed

$(SHLIB): libmseed/libmseed.a

libmseed/libmseed.a:
        @(cd libmseed && $(MAKE) static CC="$(CC)" CFLAGS="$(CFLAGS)")

[03:20]

4) Convert example C code to R callable routine.

This step assumes that you already know everything learned in the .Call(“hello”) post.  I took the code in msview.c and whittled it down as much as possible before converting it into an R callable function that spits out, once again, “Hello World!”.  Here’s what I ended up with.  Copy this code to mseed/src/test_hello_mseed.c

/***************************************************************************
 R callable function
 ***************************************************************************/

#include <R.h>
#include <Rdefines.h>

#include <libmseed.h>

SEXP test_hello_mseed () {

  SEXP result;

  static flag  verbose    = 0;
  static int   reclen     = -1;

  MSRecord *msr = 0;

  int dataflag   = 0;
  int64_t totalrecs  = 0;
  int64_t totalsamps = 0;
  int retcode;

  char *inputfile = "mseed/src/libmseed/example/test.mseed";

  /* Loop over the input file */
  while ( (retcode = ms_readmsr (&msr, inputfile, reclen, NULL, NULL, 1,
                 dataflag, verbose)) == MS_NOERROR )
    {
      totalrecs++;
      totalsamps += msr->samplecnt;

      printf ("Hello World!\n");
    }

  if ( retcode != MS_ENDOFFILE )
    ms_log (2, "Cannot read %s: %s\n", inputfile, ms_errorstr(retcode));

  /* Make sure everything is cleaned up */
  ms_readmsr (&msr, NULL, 0, NULL, NULL, 0, 0, 0);

  PROTECT(result = NEW_INTEGER(2));
  INTEGER(result)[0] = (int) totalrecs;
  INTEGER(result)[1] = (int) totalsamps;
  UNPROTECT(1);

  return(result);
}

[04:30]

You may have noticed that we hardcoded the location of our data file.  Not really following best practices here but we’re just trying to do the minimum to get it all to work.  We just need to make sure we’re in the directory above mseed/ when we invoke this function.

5) Add wrapper code, NAMESPACE and DESCRIPTION

As described in  .Call(“hello”) we need to add a wrapper function for our C code.  Create mseed/R/mseedWrappers.R with the following content:

# wrapper function to invoke test_hello_mseed
test_hello_mseed <- function() {
  result <- .Call("test_hello_mseed")
  return(result)
}

Now we need to educate R about the namespace associated with our package.  Create mseed/NAMESPACE with the following:

# Import required packages
# none so far

# Load dynamic libraries (shared object files)
useDynLib(mseed)

# Export functions defined in R code
export("test_hello_mseed")

When R compiles everything, the test_hello_mseed() function written in C will be part of the mseed shared object library.

The last step needed for package compilation is the mseed/DESCRIPTION file:

Package: mseed
Type: Package
Title: Test wrapper for libseed
Version: 0.0-0
Date: 2012-12-12
Author: Jonathan Callahan <jonathan.s.callahan@gmail.com>
Maintainer: Jonathan Callahan <jonathan.s.callahan@gmail.com>
Description: A package that demonstrates how to wrap a C library using
             the seismic data Mini-SEED library as an example
            (http://www.iris.edu/software/libraries/).
License: GPL (>= 2)
Depends: R (>= 2.14)
Collate: mseedWrappers.R

[06:50]

6) Wave your wand and watch the magic!

Assuming you have everything in place and have copy-pasted well the following should work if you position yourself above the mseed/ directory.  (Always start with R CMD REMOVE mseed to get rid of any previous version.)

$ R CMD REMOVE mseed
Removing from library ‘/usr/lib/R/library’
Updating HTML index of packages in '.Library'
Making packages.html  ... done
$ R CMD build mseed
* checking for file ‘mseed/DESCRIPTION’ ... OK
* preparing ‘mseed’:
* checking DESCRIPTION meta-information ... OK
* cleaning src
* checking for LF line-endings in source and make files
* checking for empty or unneeded directories
* building ‘mseed_0.0-0.tar.gz’

$ R CMD INSTALL mseed
* installing to library ‘/usr/lib/R/library’
* installing *source* package ‘mseed’ ...
** libs
gcc -m32 -std=gnu99 -I/usr/include/R -Ilibmseed -I/usr/local/include    -fpic  -O2 -g -pipe -Wall 
-Wp,-D_FORTIFY_SOURCE=2 -fexceptions -fstack-protector --param=ssp-buffer-size=4 -m32 -march=i686 
-mtune=atom -fasynchronous-unwind-tables -c test_hello_mseed.c -o test_hello_mseed.o
make[1]: Entering directory `/home/jonathan/mseed/src/libmseed'
gcc -m32 -std=gnu99 -O2 -g -pipe -Wall -Wp,-D_FORTIFY_SOURCE=2 -fexceptions -fstack-protector 
--param=ssp-buffer-size=4 -m32 -march=i686 -mtune=atom -fasynchronous-unwind-tables   
-c -o fileutils.o fileutils.c
...
** building package indices ...
** testing if installed package can be loaded

* DONE (mseed)
Making packages.html  ... done
$ R --vanilla
...
> library(mseed)
> test_hello_mseed()
Hello World!
...
Hello World!
[1]   36 4200
>

[08:20]

It’s a miracle!!!

Congratulations are in order.  Go have a coffee or a beer or whatever to kill the rest of your 15 minutes … unless …

Unless, that is, you want to tackle actually communicating data between R and this C library.  That will take a few more minutes

7) Start with a C code example that reads data from a buffer

Presumably, we’re packaging our library so that R can pass objects to C code and get objects back.  That means that our C code better know how to read data from a buffer.  We probably want to start out with a C code example that does this.  In the real world you will either have to write the buffer example yourself or ask a friend, preferably the author of the package, to help out.  I started with msviewmemory.c.  Better make sure it works:

$ cd mseed/src/libmseed/example
$ wget http://mazamascience.com/Published/msviewmemory.c
--2012-11-16 16:42:11--  http://mazamascience.com/Published/msviewmemory.c
...
$ make msviewmemory
cc -I..  -L..  msviewmemory.c  -lmseed -o msviewmemory
$ ./msviewmemory test.mseed
IU_COLA_00_LHZ, 000001, M, 512, 112 samples, 1 Hz, 2010,058,06:50:00.069539
...

[09:30]

So far so good.

8) Whittle this code down to an R callable routine

You can try and figure out what the code below does later.  We’re at risk of going over time so go back to the directory above mseed and just copy and paste the following code into mseed/src/test_hello_mseed_mem.c:

/*************************************************************************** 
 R callable function to pass raw data to a C test function
 ***************************************************************************/
 #include <R.h>
#include <Rdefines.h>

#include <stdio.h>
#include <sys/stat.h>
#include <libmseed.h>

SEXP test_hello_mseed_mem (SEXP buffer) {

  int bufferLength;
  char *bufferPtr;
  SEXP result;

  PROTECT(buffer = AS_RAW(buffer));
  bufferPtr = RAW_POINTER(buffer);
  bufferLength = LENGTH(buffer);

  for (int i=0; i<10; i++) {
    Rprintf("buf[%d]=0x%02x\n",i,bufferPtr[i]);
  }

  Rprintf("bufferLength = %d\n",bufferLength);

  long long int bsize = (long long int) bufferLength;
  long long int boffset = 0;

  static flag  verbose    = 0;
  static int   reclen     = -1;

  MSRecord *msr = 0;

  int dataflag   = 0;
  int64_t totalrecs  = 0;
  int64_t totalsamps = 0;

  /* Loop over the buffer */
  boffset = 0;
  while ( boffset < bsize )
    {
      if ( msr_parse (bufferPtr+boffset, bsize-boffset, &msr, reclen, dataflag, verbose) )
    {
      if ( verbose )
        ms_log (2, "Error parsing record at offset %lld\n", boffset);

      boffset += 256;
    }
      else
    {
      totalrecs++;
      totalsamps += msr->samplecnt;

          printf ("Hello World!\n");

          boffset += msr->reclen;
    }
    }

  /* Make sure everything is cleaned up */
  ms_readmsr (&msr, NULL, 0, NULL, NULL, 0, 0, 0);

  PROTECT(result = NEW_INTEGER(3));
  INTEGER(result)[0] = (int) bufferLength;
  INTEGER(result)[1] = (int) totalrecs;
  INTEGER(result)[2] = (int) totalsamps;

  UNPROTECT(2);
  return(result);
}

[11:10]

9) Add a wrapper, external data and a demo

Add this to mseed/R/mseedWrappers.R:

# wrapper function to invoke test_hello_mseed_mem
test_hello_mseed_mem <- function(buffer) {
  result <- .Call("test_hello_mseed_mem",buffer)
  return(result)
}

Of course, we’ll also have to add one more line specifying this new function to mseed/NAMESPACE:

export("test_hello_mseed_mem")

Make a copy of test.mseed so that it is found in the normal place for package “external data”:

$ cp mseed/src/libmseed/example/test.mseed ./mseed/inst/extdata/

Now add the following demonstration script as mseed/demo/hello.R:

# Demonstrate successful processing of a miniseed file

library(mseed)

# Read in a mseed file as raw bytes
mseedFile <- system.file("extdata", "test.mseed", package="mseed")
rawData <- readBin(mseedFile, raw(), n=1e+6)

# How many bites do we have?
length(rawData)

# Let's look at the first 10
rawData[1:10]

# Now pass the data to test_hello_mseed_mem()
test_hello_mseed_mem(rawData)

and a file named mseed/demo/00Index:

hello           read mseed file and pass it to a C test function

[12:55]

10) Racing to the finish line

$ R CMD REMOVE mseed
...
$ R CMD build mseed
...
$ R CMD INSTALL mseed_0.0-0.tar.gz 
...
]$ R --vanilla
...
> library(mseed)
> demo(package="mseed")
> demo(hello)
...
Hello World!
[1] 18432    36  4200
>

[Smack that timer!] Done.

So, how is your time?  Did you keep up the pace and finish in under 15 minutes?  Or did you get distracted trying to understand what you were doing and go over time.

You should feel good about your accomplishment in either case if you got this far.  (I certainly did!)  Compiling a C library for inclusion with R is not for the faint of heart.  I hope that this baby-steps example makes it seem like it’s at least possible if perhaps not easy.

Best of luck extending R!

In two previous posts we described how R can call C code with .C() and the more complex yet more robust option of calling C code with .Call().  Here we will describe how the Rcpp package can be used to greatly simplify your C code without forcing you to become expert in C++.

First off, kudos to Dirk Eddelbuettel and Romain François who’s tireless efforts to improve, promote and document Rcpp have produced one of CRAN’s most popular packages.  As of Rcpp version 0.9.15 there are 82 “Reverse depends” — other packages that utilize Rcpp.  There are also eight vignettes that describe the package in human terms.  Even more accessible are Dirk’s papers and presentations and Romain’s blog.  While we’re heaping praise, let’s not forget to mention Romain’s graph gallery.

Even after all that praise and all the available documentation we’re still left with the problem of where to start.  There are no vignettes targeted at the R user who wants to try out a couple of C routines but isn’t otherwise inclined to learn C++ and who doesn’t want to write an entire package — at least not yet.  A good place to review the motivation behind Rcpp and see some great example code would be Dirk’s 62 slide  Rcpp: Seamless R and C++ Integration presentation or the longer, 146 slide Rcpp Tutorial Parts I and II.  Between those two presentations and the package documentation you really have access to all the information you need.  However, in the interest of continuing our “baby steps” theme, we will once again recreate our three “Hello World!” examples, this time using Rcpp.

Getting Set Up

The very first thing you will have to do is install the Rcpp package.  (See Installing Packages if this is unfamiliar.)  If you keep your version of R at the bleeding edge, you can do this by invoking R and telling it to just grab the latest package.

$ R --vanilla

R version 2.14.1 (2011-12-22)
Copyright (C) 2011 The R Foundation for Statistical Computing
...
> install.packages( repos=c('http://cran.fhcrc.org/'), pkgs=c('Rcpp') )
Installing package(s) into ‘/usr/lib/R/library’
(as ‘lib’ is unspecified)
Warning message:
In getDependencies(pkgs, dependencies, available, lib) :
  package ‘Rcpp’ is not available (for R version 2.14.1)

Oops!  Some of us are slightly behind the times.  Rather than upgrade my version of R today, I’ll instead take a look inside the Rcpp Old sources: archive and find a version that was released close to the (2011-12-22) date of my version of R.  Version 0.9.10 was released on 17-Feb-2012 and should be compatible.  We’ll install that version from the command line with the following:

$ wget http://cran.r-project.org/src/contrib/Archive/Rcpp/Rcpp_0.9.10.tar.gz
...
$ sudo R CMD INSTALL Rcpp_0.9.10.tar.gz
[sudo] password for jonathan: 
* installing to library ‘/usr/lib/R/library’
* installing *source* package ‘Rcpp’ ...
...
** testing if installed package can be loaded

* DONE (Rcpp)
Making packages.html  ... done

The last thing we need to do is let our compilers and linkers know where the new, Rcpp libraries are located (see slide 48 of the Tutorial).

$ export PKG_CPPFLAGS=`Rscript -e "Rcpp:::CxxFlags()"`
$ export PKG_CPPFLAGS=`Rscript -e "Rcpp:::LdFlags()"`

Yes, Whew! once again.  But now we are ready to write C code that looks more like C code and, once we embrace Rcpp’s syntactic sugar, perhaps even like R code.

Baby steps example

When using the Rcpp package we are still using the .Call() interface to C code.  All of the changes will be seen in the C code which now becomes C++ code with a .cpp extension.  Because C++ is a superset of C, this is perfectly legal but we will need to educate Rcpp by using RcppExport.  According to slide 136 of the Tutorial:

* note : RcppExport is an alias to ‘extern "C"‘ defined by Rcpp.
*
* It gives C calling convention to the [...] function so that
* it can be called from .Call in R. Otherwise, the C++ compiler mangles the
* name of the function and .Call can’t find it.

OK.  Here we go with our first  C++ example: helloA2.cpp.

#include <Rcpp.h>
RcppExport SEXP helloA2() {
  printf("Hello World!\n");
  return(R_NilValue);
}

It looks remarkably similar to helloA1.c from the previous post and is compiled and invoked in much the same way.

R CMD SHLIB helloA2.cpp
g++ -m32 -I/usr/include/R -L/usr/lib/R/library/Rcpp/lib -lRcpp -Wl,-rpath,/usr/lib/R/library/Rcpp/lib 
-I/usr/local/include   -I/usr/lib/R/library/Rcpp/include -fpic  -O2 -g -pipe -Wall 
-Wp,-D_FORTIFY_SOURCE=2 -fexceptions -fstack-protector --param=ssp-buffer-size=4 -m32 -march=i686 
-mtune=atom -fasynchronous-unwind-tables -c helloA2.cpp -o helloA2.o
g++ -m32 -shared -L/usr/local/lib -o helloA2.so helloA2.o -L/usr/lib/R/library/Rcpp/lib -lRcpp 
-Wl,-rpath,/usr/lib/R/library/Rcpp/lib -L/usr/lib/R/lib -lR

The R wrapper code is identical to the one in the previous post for helloA1.c as is the usage in an R session.

# call RCpp C code
dyn.load("helloA2.so")
helloA2 <- function() {
  result <- .Call("helloA2")
}

> source('wrappers.R')
> greeting <- helloA2()
Hello World!
> class(greeting)
[1] "NULL"

Simpler C code

We’ll leave out the wrappers and R session in the next two examples as they are identical to the examples in the previous post.  But just look at how much simpler the C code gets!  Instead of allocating memory, protecting from garbage collection and all that casting between types we get this for helloB2.cpp:

#include <Rcpp.h>
RcppExport SEXP helloB2() {
  Rcpp::StringVector result(1);
  result[0] = "Hello World!";
  return(result);
}

and this for helloC2.cpp:

#include <Rcpp.h>
RcppExport SEXP helloC2(SEXP greetingPointer) {
  Rcpp::StringVector greeting(greetingPointer);
  Rcpp::NumericVector result(greeting.size());
  for (int i=0; i<greeting.size(); i++) {
    result[i] = strlen(greeting[i]);
  }
  return(result);
}

Now we’re talkin’.  This is starting to look like code a C programmer, Heck, even a Java programmer could get comfortable with.  The tools provided by Rcpp are systematic, readable and, as a welcome change, very well documented.  Whereas I was hesitant to recommend writing C code for the .Call() interface because of the painful learning curve, I am happy to report that the learning curve for using Rcpp does not require that you first climb a mountain to learn C++.  C programmers of all skill levels will benefit from using Rcpp.

More Examples

While a “Hello World!” example may be a great place to start, it is unlikely to provide any useful template code for people.  For that you will want to poke around in the source code of the many packages that depend on Rcpp.  Dirk Eddelbuettel is of course quite interested in these dependent packages and describes some of them on slide 42 of Seamless R and C++ Integration and slide 29 of Rcpp Tutorial Parts I and II.  I was pleased to learn that some of these packages are written in C and will hopefully provide excellent example code.

Hadley Wickham has written a comprehensive tutorial for the Rcpp package.

To everyone who is trying to improve and extend R — Best of Luck!