Using R: Working with Geospatial Data (and ggplot2)

Loading libraries and data

# load libraries
library(ggplot2)
library(sp)
library(rgdal)
library(rgeos)

# create a local directory for the data
localDir <- "R_GIS_data"
if (!file.exists(localDir)) {
  dir.create(localDir)
}

# download and unzip the data
url <- "ftp://www.ecy.wa.gov/gis_a/inlandWaters/wria.zip"
file <- paste(localDir, basename(url), sep='/')
if (!file.exists(file)) {
  download.file(url, file)
  unzip(file,exdir=localDir)
}

# create a layer name for the shapefiles (text before file extension)
layerName <- "WRIA_poly"

# read data into a SpatialPolygonsDataFrame object
dataProjected <- readOGR(dsn=localDir, layer=layerName)

Transforming the data

Thus far, we haven’t done anything radically different than before, but in order to prepare the data for plotting in a ggplot, we’ll have to do a couple manipulations to the structure of the data. ggplot2 will only work with a data.frame object, so our object of class of SpatialPolygonsDataFrame will not be appropriate for plotting. Let’s write some code and discuss why this kind of transformation is necessary.

# add to data a new column termed "id" composed of the rownames of data
dataProjected@data$id <- rownames(dataProjected@data)

# create a data.frame from our spatial object
watershedPoints <- fortify(dataProjected, region = "id")

# merge the "fortified" data with the data from our spatial object
watershedDF <- merge(watershedPoints, dataProjected@data, by = "id")

# NOTE : If we so choose, we could have loaded the plyr library to use the
#      : join() function. For those familiar with SQL, this may be a more
#      : intuitive way to understand the merging of two data.frames. An
#      : equivalent SQL statement might look something like this:
#      : SELECT *
#      : FROM dataProjected@data
#      : INNER JOIN watershedPoints
#      : ON dataProjected@data$id = watershedPoints$id

# library(plyr)
# watershedDF <- join(watershedPoints, dataProjected@data, by = "id")

What does all this code mean and why do we need it? Let’s go through this line by line.

dataProjected@data$id <- rownames(dataProjected@data)

Here we are appending to the data an extra column called “id”. This column will contain the rownames so that we define an explicit relationship between the data and the polygons associated with that data.

watershedPoints <- fortify(dataProjected, region = "id")

Fortify? What does that even mean? A quick search on the internet will yield some helpful documentation. (See fortify.sp documentation). Basically, fortify take two arguments: model, which will consist of the SpatialPolygonsDataFrame object we wish to convert and region, the name of the variable by which to split regions. If all goes according to plan, some magic happens and we get a data.frame, just like we wanted… well, not quite. If you inspect this data.frame, you’ll notice it appears to be missing some critical information. Fret not! Using the relationship we created earlier, we can merge these two datasets with the following command.

watershedDF <- merge(watershedPoints, dataProjected@data, by = "id")

And viola! Now that we’ve created a data.frame that ggplot2 likes, we can begin plotting. Before we get to plotting, let’s take a quick look at this new data.frame we’ve created.

head(watershedDF)