Preface

Aim

Have you ever wanted to create a fictitious map? Perhaps you wanted to test a new spatial model on simulated data, and didn’t want to use real spatial units. The purpose of this document is to demonstrate how to create fictitious maps that either resemble artificial square grids or more realistic looking geographic units like counties/suburbs. Code for exporting these maps as ESRI shapefiles is also provided.

I plan to update this guide in the future. So if you find any mistakes, have any suggestions for amendments or additions, or have any questions about the content, please contact me at earl.w.duncan@gmail.com.

R Libraries

The R code in this document requires the following packages to be installed and loaded:

library(sp)             # For coordinates(); Polygons(); SpatialPolygons(); SpatialPolygonsDataFrame()
library(ggplot2)        # For fortify(); ggplot()
library(rgdal)          # For writeOGR()

1 Regular Lattice (Square Grid Graph)

1.1 Understanding the Components

A square grid graph is a simple lattice graph consisting of square polygons. It’s simplicity can be attractive as it is easy to construct, fast to render, readily conveys fundamental concepts in topological studies such as autocorrelation and spatial patterns, and doesn’t detract attention from other more important aspects of an application such as the values/colours that may be mapped to these squares.

To understand what data is required in the construction of a square grid graph, let’s begin by examining the first 7 lines of a data.frame representation of a square grid graph.

##   long lat order  hole piece id group
## 1    0   0     1 FALSE     1  1   1.1
## 2    0   1     2 FALSE     1  1   1.1
## 3    1   1     3 FALSE     1  1   1.1
## 4    1   0     4 FALSE     1  1   1.1
## 5    0   0     5 FALSE     1  1   1.1
## 6    1   0     1 FALSE     1  2   2.1
## 7    1   1     2 FALSE     1  2   2.1

The first 5 rows relate to the first square (bottom-left of the lattice). The first two columns give the coordinates of each vertex (note that the first coordinate, 0, 0 is repeated as the final coordinate to close the polygon). The next column, order, specifies the order these vertices are joined to create the polygon. Note that in the popular ESRI shapefile data format, it is important that polygon vertices are connected in a clockwise direction, unless the polygons represent holes in which case they are connected in an counter-clockwise direction. This brings us to the fourth column, hole, which specifies whether the polygon represents a hole or not (the same value, either TRUE or FALSE, would be the same for all 5 rows relating to the given polygon). The next column, piece is used to identify multiple polygons which should be treated as a single Polygons entity. I explain this in a bit more detail in another guide. For the problem at hand, this variable isn’t important – this value will always be 1 for each square. id is simply a unique identifier for each square – all vertices relating to the same square will have the same id. Lastly, group is just a concatenation of the id and piece variables to facilitate discriminating between individual polygons. Since piece will always take the value 1, this is synonymous to the id variable.

1.2 Construction

To construct a square grid graph, we need to create a data.frame representation of the graph like the one shown above, and then convert it to a SpatialPolygonsDataFrame object. This is achieved by the following user-defined function (UDF) which takes the user-specified inputs nrow and ncol which denote the number of rows (height) and columns (width) of the lattice respectively.

# UDF which creates a regular grid of size nrow by ncol
make.lattice <- function(nrow, ncol){
    if(missing(ncol)){
        ncol <- nrow
    }
    # Create coordinate data
    id <- rep(1:(nrow * ncol), each = 5)
    x <- c(0, 1, 1, 0, 0) + (id-1) %% ncol
    y <- c(0, 0, 1, 1, 0) + rep((1:nrow)-1, each = ncol*5)
    coords <- data.frame(id, x, y, order = 1:5, hole = FALSE, piece = 1, group = id)
    
    # Convert coordinates from data.frame to SpatialPointsDataFrame
    coordinates(coords) <- ~x + y
    
    # Convert SpatialPointsDataFrame to a SpatialPolygonDataFrame
    ID <- data.frame(id = unique(id))
    lattice <- points2polygons(coords, ID)
    return(lattice)
}

If ncol isn’t specified, then the lattice defaults to a square (of squares), with nrow rows and columns.

The four lines leading up to and including coords construct a data.frame like the one shown above. At this stage, the set of points defined by each row don’t actually form polygons. To do this, we first convert these points to a SpatialPointsDataFrame using the coordinates function from the sp package. This is subsequently converted to a SpatialPolygonDataFrame using the UDF points2polygons, defined below.

# UDF for converting spatial points to spatial polygons
points2polygons <- function(spdf, data, df.out = TRUE){
    get.grpPoly <- function(group, ID, spdf){
        Polygon(coordinates(spdf[spdf$id == ID & spdf$group == group,]))
    }
    get.spPoly  <- function(ID, spdf){
        Polygons(lapply(unique(spdf[spdf$id == ID,]$group), get.grpPoly, ID, spdf), ID)
    }
    spPolygons  <- SpatialPolygons(lapply(unique(spdf$id),get.spPoly, spdf))
    if(df.out){
        row.names(data) <- data$id
        spPolygons <- SpatialPolygonsDataFrame(spPolygons, match.ID = TRUE, data = data)
    }
    return(spPolygons)
}

This UDF is quite intricate, and the specifics aren’t important for the purpose of this guide. It suffices to say that this UDF performs the conversion we desire.

1.3 Example

Once the UDFs above are defined, creating a regular lattice is very straightforward – simply pass the desired input values to the make.lattice function.

# Choose grid size (N = nrows, M = ncols)
N <- 4
M <- 8

# Create the lattice
map <- make.lattice(N, M)

# Plot the lattice using base R graphics
plot(map)

# Plot the lattice using ggplot2
map.df <- fortify(map)
ggplot(map.df, aes(x = long, y = lat, fill = group)) +
    geom_polygon() +
    coord_map() +
    scale_fill_discrete(guide = FALSE)

I deliberately named the resulting object map because this is how such objects are typically used in practice, and it lends itself to the obvious extension of forming maps out of irregular polygons instead of a regular grid.

1.4 Exporting

To save this spatial object to disk, use writeOGR. The following example demonstrates how to save it in the commonly used ESRI shapefile format.

# Not run
fp.wd <- getwd()    # Or specify the desired directory
writeOGR(map, dsn = fp.wd, layer = "My Square Lattice", driver = "ESRI Shapefile")

2 Irregular Polygons

This is a work in progress. I plan to update this at a later date after more research and testing.

Creating Artificial Maps and Shapefiles in R