Getting rasters into shape from R

Today I needed to convert a raster to a polygon shapefile for further processing and plotting in R. I like to keep my code together so I can easily keep track of what I’ve done, so it made sense to do the conversion in R as well. The fantastic raster package has a function that can do this (rasterToPolygons), but it seems to take a very, very long time and chew up all system resources (at least on my not particularly souped up X200 laptop).

We can cut down the time taken for conversion by calling a GDAL utility,  gdal_polygonize.py, directly from R using system(). GDAL needs to be installed on your system, but you’ll probably want it installed anyway if you’re planning to talk to spatial data (in fact, you might find you actually already have it installed, as it is required for a bunch of other GIS software packages). Together with it’s cousin OGR, GDAL allows reading and writing of a wide range of raster and vector geospatial data formats. Once installed, you have at your disposal a bunch of GDAL Utilities that can be run from the terminal, including gdal_polygonize.py. For polygonize and a couple of others, you’ll also need Python installed (again, you may already have it installed, so check first!).

Anywho… let’s compare rasterToPolygons with gdal_polygonize. The function included below (gdal_polygonizeR) borrows from code provided in a post by Lyndon Estes on R-sig-geo. Thanks mate!

R functions to filter rjags results

A while back I was running a bunch of JAGS models through R, using the rjags (written by Martyn Plummer) and R2jags (by Yu-Sung Su) packages. These packages provide a great interface to the JAGS software, which allows analysis of Bayesian models (written in the BUGS language) through Markov chain Monte Carlo simulation.

Running a JAGS model using these tools returns an rjags object, which when printed to the screen, summarises the posterior distribution of each monitored node, giving its mean and standard deviation, a range of quantiles, and its Gelman-Rubin convergence diagnostic statistic (Rhat), which indicates the ratio of variance within chains to that among chains. The summary is great, but when monitoring a large number of nodes, printing these to the screen can cause R to hang, and can exceed the screen buffer (not to mention making it painful to find the nodes you’re immediately interested in).

To help deal with this I wrote a couple of simple R functions:

• jagsresults: return a matrix containing summary results for just the nodes you are interested in (using regular expression pattern-matching, if desired).
• rhats: sort the output by the nodes’ Rhat values, making it easy to show the n least converged nodes.

I, Rbot: Tweeting from R

Model status update. A Tweet from the Rbot.

Over the past few weeks I’ve been running batches of JAGS simulations from R. Although these models typically converged within an hour or so, more complex models can take days, or even weeks to converge. Because we, as humans, are required to bathe, feed, socialise and attend weddings, there will inevitably come a time when you are required to leave the safety of your modelling den while your simulation is still chugging away. It would be nice, however, to keep in touch with your models’ progress, to be made instantly aware of any errors that occur in your absence, and to be advised upon successful task completion. This would be decidedly more satisfying than returning from a week-long holiday to find that your model broke on Day 1, just after you closed the door on your way out.

Enter Twidge, a command line Twitter client that humbly fills this neglected niche by allowing you to send a short Tweet to yourself from R via the system() function. The tweet can, of course, be a message pasted together from R objects, which permits dynamic tweet content and means your R-Tweeting power is really only limited by your imagination (and the 140 character message cap).

OpenBUGS error messages. Something went right.

Although my relationship with OpenBUGS is still in the early stages, one glaringly obvious and quite infuriating shortcoming has been the lack of detail in its error messages. The initiated amongst you are no doubt familiar with the phrase “Sorry, something went wrong in procedure [blah]”. These error messages had me chasing my tail for hours in vain attempts to find offending segments of code. Surely, I thought, there must be a way to display more detailed information about the cause of an error.

A default OpenBUGS error message

As it happens, there is.