Executive Summary

The following code will help you build your own maps in R using base plotting, Lattice plot methods for spatial data, the ggplot2 system, the GoogleVis Chart API and interactive javascript visualizations. Enjoy and have fun with it!


You may want to download the following files:

  1. First, you need a shapefile of the country you want to build a map for. The shapefile format is a popular geospatial vector data format for geographic information system (GIS) software. It is developed and regulated by Esri. For our examples, we will use two shapefiles for Norway: The first shapefile of Norway is available from the website of the Norwegian Kartverket. You can download the file here. We will use the shapefile called NO_Arealdekke_pol.shp. Don’t forget to read the terms of use.

  2. The second shapefile is available from the website of the Global Administrative Areas at this location. Download the shapefile format for Norway. You can also download shapefiles for other countries there.

  3. You may also need a file containing the X-Y coordinate of spatial locations you want to plot. In our case, we will use a file containing the Norwegian postal code areas. This file is available here. I have also used the coordinates of cities and other areas from this website.

Building the Map of Norway in R

First we load the necessary R packages, download the needed files and set up the directory.

# Installing all the libraries to build maps in R
library(plyr)         # To manipulate data
library(ggplot2)      # To have ggplot2 graphic interface
library(lattice)      # To have Lattice graphic interface
library(rgdal)        # To load "shapefiles" into R and use in conversions of spatial formats 
library(rgeos)        # To use in geometrical operations
library(spatstat)     # To use with Voronoi Tesselation
library(sp)           # Methods for retrieving coordinates from spatial objects.
library(maptools)     # A package for building maps
library(maps)         # A package for building maps
library(RColorBrewer) # To build RColorBrewer color palettes
library(grDevices)    # To build grDevices color palettes
library(reshape2)     # To have more flexibility when transforming data
library(rCharts)      # To create and customize interactive javascript visualizations in R
library(knitr)        # For dynamic report generation in R
library(base64enc)    # Tools for base64 encoding
suppressPackageStartupMessages(library(googleVis)) # To use with Google visualization in R

# Setting your working directory
setwd("/Users/isa/Documents/Courses/Johns Hopkins/RKart_Norge")

Building Maps using R’s Base Plotting System

There are several ways to build maps in R. First let’s build a simple map of Norway showing a few Norwegian cities. We will use R’s base plotting system.

# First, let's load the first shapefile of Norway into R
fylke <- readOGR(dsn="/Users/isa/Documents/Courses/Johns Hopkins/RKart_Norge/33_N5000_shape" ,
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/isa/Documents/Courses/Johns Hopkins/RKart_Norge/33_N5000_shape", layer: "NO_Arealdekke_pol"
## with 1324 features
## It has 5 fields

# Secondly, we create a simple dataframe with a few Norwegian cities coordinates
name <- c("Postnummer", "LAT", "LON"); Oslo <- c("0252", 59.9084400, 10.7230100); 
Larvik <- c("3252",59.0686800,10.0473000); Bergen <- c("5052",60.3713600,5.3537800);  
Haugesund<- c("5532",59.4094000,5.2995100); Trondheim <- c("7052",63.4240100,10.4409600);
Nordkapp <- c("9751",70.9821400,25.9707600); Tromso <- c("9252",69.6478700,18.9566000)
cities <- as.data.frame(rbind(Oslo, Larvik, Bergen, Haugesund, Trondheim, Tromso, Nordkapp),col.names=name, stringsAsFactors=FALSE) 
cities$V1 <- as.character(cities$V1); cities$V2 <- as.numeric(cities$V2); cities$V3 <- as.numeric(cities$V3)
colnames(cities) <- name;

# We check classes and formats
proj=CRS(proj4string(norway)); proj; class(cities) 
## CRS arguments:
##  +proj=utm +zone=33 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0
## [1] "data.frame"
# We convert the dataframe "cities" to a spatial object
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"
# The spTransform function provides transformation between datum(s) and projections
cities.proj=spTransform(cities, proj)
# We check that norway and cities.proj are identical before merging them
## [1] TRUE
# We are now ready to plot the first map with R'base plotting system.
plot(norway, border="darkgrey")
# We plot the cities coordinates
points(cities.proj, col="red", pch=16)
# We want to add the city names
text(cities.proj, labels=row.names(cities), cex=0.6, font=2, offset=0.5, adj=c(0,2))
# We need a legend
legend("bottomright", legend=c("©Kartverket"),col=black, cex=0.5, bty="n")

Performing a Voronoi Tesselation with Base Plotting in R

Now, let’s build a map showing all postal code areas with a Voronoi Tesselation:

  1. We plot the map of Norway as previously shown with the first shapefile of Norway.

  2. We plot all the 4463 postal codes coordinates just as we plotted a few cities earlier.

  3. We perform a Voronoi Tesselation. See the following wiki link for more details on Voronoi Tesselation. We use the “spatstat” package in R to do this.

# We load the postal codes into R.
postal <- read.delim("postnummer.csv", colClasses="character")
postal <- transform(postal, LON=as.numeric(LON), LAT=as.numeric(LAT))
# We remove all postal codes not in use
postal2=postal[!grepl("ikkje i bruk", postal$POSTSTAD),]
# Some locations have the same postal codes. We transform the data to avoid plotting problems. 
post=ddply(postal2, ~LON+LAT, summarise, POSTNR=paste(POSTNR, collapse=", "))
# We convert to the right spatial format
proj=CRS(proj4string(norway)); coordinates(post)=~LON+LAT
proj4string(post)=CRS("+init=epsg:4326");  post.proj=spTransform(post, proj)
# We create a map of Norway and plot the postal codes
plot(norway, border="darkgrey")
points(post.proj, col="red", pch=".")

We are now ready to perform the Voronoi Tesselation.

# We perform a Voronoi Tesselation using the "spatstat" package in R.
bboks.owin=owin(bboks.norway[1,], bboks.norway[2,])
post.ppp=as.ppp(coordinates(post.proj), W=bboks.owin)
diagram.poly=as(diagram, "SpatialPolygons")
 # We plot the results with base plotting