Plotting spatial blocks borders

plot_blocks.R

# This procedure allows to plot blocks borders
# over some spatially-arranged variable

# Data is assumed to be on a regular grid
# The dataset contains row and column coordinates
# a numerical variable indicating the block number of the current cell
# and the numerical variable to be plotted

# Requirements:
# either package rgeos, which depends on [GEOS](http://trac.osgeo.org/geos/) (GEometry Open Source)
# either package gpclib (which interfaces some non-foss libraries)
# Recommended for Windows/Mac users is gpclib, *but* you need to install from **source**
#   i.e. need [RTools](http://cran.r-project.org/bin/windows/Rtools/) to compile some C code

install.packages(gpclib, type = 'source')

library(sp)
library(maptools)

dat <- data.frame(expand.grid(1:9, 1:9),
                  c(rep(1:3, 3, each=3),
                    rep(4:6, 3, each = 3),
                    rep(7:9, 3, each = 3)),
                  runif(81))
names(dat) <- c('x', 'y', 'bl', 'z')

matrix(dat[order(dat$x, dat$y), 'bl'], 9, 9)

dat.grid <- dat
coordinates(dat.grid) <- 1:2   # promote to SpatialPointsDataFrame
gridded(dat.grid) <- TRUE      # promote to SpatialPixelsDataFrame

image(dat.grid['z'])

poly <- as(dat.grid, 'SpatialPolygons')

if (!rgeosStatus()) gpclibPermit()
blpoly <- unionSpatialPolygons(poly, IDs = dat$bl)

plot(blpoly, add = TRUE, border ='grey', lwd = 2)