is a free software environment and language for statistical analysis and data visualisation. Its programming language is open source and many statisticians and data scientists write packages for the software. R usage is increasing exponentially and many major companies like Google and Facebook run complex statistical analyses with it.

• Open source: No cost
• Performance: R is stable and blazingly fast
• Reproducibility: code can be shared and the results repeated
• Transparency: code can be checked forensically and any errors spotted
• Automation: statistical analyses and data visualisations can be run and re-run with new and existing data sets
• Support network:
  • books: e.g. Brunsdon & Comber (2015)
  • online tutorials: e.g. spatial.ly; Creating-maps-in-R
  • online forums: e.g. StackOverflow; GIS.StackExchange
  • and cheatsheets
• Integrated with ArcGIS: ArcGIS Bridge is an R package which allows ArcGIS and R to dynamically access data
• and you can produce some amazing maps…….

• rgdal - reading and writing spatial data
• sp - classes and methods for spatial data
• ggplot2 - graphics package
• ggmap - extends ggplot2 with raster layers e.g. Google Maps and Stamen Maps
• leaflet - interactive web maps with the JavaScript 'Leaflet' library
• classInt - class intervals
• RColorBrewer - colour schemes for maps

For more information about spatial data packages in R visit: http://cran.r-project.org/web/views/Spatial.html

Read the London borough boundary shapefile downloaded from the London Datastore using the readOGR function from the rgdal package and create a spatial object called 'boroughs'.

boroughs <- readOGR(".", "London_Borough_Excluding_MHW", verbose = FALSE)

Inspect the attribute table

head(boroughs@data, 4)

##                   NAME  GSS_CODE  HECTARES NONLD_AREA ONS_INNER
## 0 Kingston upon Thames E09000021  3726.117      0.000         F
## 1              Croydon E09000008  8649.441      0.000         F
## 2              Bromley E09000006 15013.487      0.000         F
## 3             Hounslow E09000018  5658.541     60.755         F

and extract a particular column using the $ symbol, e.g. boroughs$NAME

Read in some publicly available police recorded crime data

df <- read.csv("crime_data.csv", header = T)
head(df, 4)

##       long      lat              category
## 1 0.137065 51.58367 Anti-social behaviour
## 2 0.140035 51.58911 Anti-social behaviour
## 3 0.140192 51.58231 Anti-social behaviour
## 4 0.134947 51.58806 Anti-social behaviour

Convert the dataframe to a SpatialPointsDataFrame called 'crimes' with latitude / longitude projection

coords <- SpatialPoints(df[,c("long","lat")])
crimes <- SpatialPointsDataFrame(coords, df)
proj4string(crimes) <- CRS("+init=epsg:4326")

Transform the borough boundary from epsg:27700 (British National Grid) to the epsg:4326 (WGS84) coordinate reference system.

boroughs <- spTransform(boroughs, CRS("+init=epsg:4326"))

Check the projection

boroughs@proj4string

## CRS arguments:
##  +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84
## +towgs84=0,0,0

Plot the London borough boundaries and the crime data

plot(crimes)
plot(boroughs, border = "blue", add = T)

Clip the crimes to the borough boundaries

proj4string(crimes) <- proj4string(boroughs)
crimes <- crimes[boroughs, ]

Plot the London borough boundaries and the crime data

plot(crimes)
plot(boroughs, border = "blue", add = T)

Join attributes from the borough boundary polygon layer to the crime data points layer.

borough_attributes <- over(crimes, boroughs[,c("NAME", "GSS_CODE")])
crimes$borough <- borough_attributes$NAME
crimes$census_code <- borough_attributes$GSS_CODE

Use the dplyr package to generate some summary statistics on Bicycle theft

crimes %>%
  as.data.frame() %>%
  filter(category == "Bicycle theft") %>%
  count(borough) %>%
  arrange(desc(n)) %>% 
  head(4)

## Source: local data frame [4 x 2]
## 
##         borough     n
##          (fctr) (int)
## 1   Westminster   143
## 2       Hackney   123
## 3 Tower Hamlets   113
## 4       Lambeth   108

Subset the Bicycle theft offences and count the number of points within each London borough using the over() function

bicycle_theft <- crimes[crimes$category == "Bicycle theft",]
proj4string(bicycle_theft) <- proj4string(boroughs)
pointsinpolygon <- over(SpatialPolygons(boroughs@polygons), 
                        SpatialPoints(bicycle_theft), returnList = TRUE)
boroughs$bicycle_theft <- unlist(lapply(pointsinpolygon, length))
head(boroughs@data, 4)

##                   NAME  GSS_CODE  HECTARES NONLD_AREA ONS_INNER
## 0 Kingston upon Thames E09000021  3726.117      0.000         F
## 1              Croydon E09000008  8649.441      0.000         F
## 2              Bromley E09000006 15013.487      0.000         F
## 3             Hounslow E09000018  5658.541     60.755         F
##   bicycle_theft
## 0            21
## 1            12
## 2            19
## 3            60

Create thematic intervals

classes <- classIntervals(boroughs$bicycle_theft, n=5, style="jenks")

Fortify the borough shapefile for plotting in ggplot2 / ggmap and merge the attribute data using the left_join() function from the dplyr package.

boroughs.f <- fortify(boroughs, region="NAME")
boroughs.f <- left_join(boroughs.f, boroughs@data, by = c("id" = "NAME"))

ggplot() +
  geom_polygon(data = boroughs.f, 
               aes(x = long, y = lat, group = group,
                   fill = bicycle_theft), color = "white", size = 0.2) +
  coord_map() +
  scale_fill_gradientn('', colours=brewer.pal(5,"RdPu"), breaks = classes$brks) +
  theme_nothing(legend = TRUE) +
  ggtitle(expression(atop("Police recorded Bicycle theft offences by borough", 
                          atop(italic("August 2015")))))

Retrieve the bounding box from the London boroughs using the bbox function from the sp package

bb <- as.vector(boroughs@bbox)
bb

## [1] -0.5103751 51.2867602  0.3340155 51.6918741

Request a Stamen Maps raster layer using the get_map function from the ggmap package

map <- get_map(location = bb, source = "stamen", maptype = "toner", crop = FALSE)

ggmap(map, extent = "device") +
  geom_polygon(data = boroughs.f, aes(x = long, y = lat, group = group,
        fill = bicycle_theft), color = "white", size = 0.2, alpha = 0.6) +
  coord_map() +
  scale_fill_gradientn('', colours=brewer.pal(5,"RdPu"), breaks = classes$brks) +
  theme_nothing(legend = TRUE) + 
  ggtitle(expression(atop("Police recorded Bicycle theft offences by borough", 
                          atop(italic("August 2015")))))

Create a popup showing the borough name and count of Bicycle theft offences

boroughs_popup <- paste0("<strong>Borough: </strong>",
                         boroughs$NAME,
                         "<br><strong>Offences (August 2015): </strong>",
                         boroughs$bicycle_theft)

Then add a colour palette for the map

colcode <- findColours(classes, c("#feebe2", "#fbb4b9", "#f768a1", 
                                  "#c51b8a", "#7a0177"))

leaflet() %>% addTiles() %>%
  addPolygons(data = boroughs,
              fillColor = colcode,
              fillOpacity = 0.6,
              color = "white",
              weight = 2,
              popup = boroughs_popup)

leaflet(boroughs) %>% addProviderTiles("CartoDB.Positron") %>% 
      fitBounds(bb[1], bb[2], bb[3], bb[4]) %>% 
      addPolygons(color = "grey", weight = 2, fillOpacity = 0.3) %>% 
      addCircleMarkers(data = bicycle_theft, ~long, ~lat, 
                       fillColor = "red", radius = 5, fillOpacity = 0.6, stroke = FALSE, 
                       clusterOptions = markerClusterOptions(
                         zoomToBoundsOnClick = TRUE,
                         spiderfyOnMaxZoom = TRUE, maxClusterRadius = 50))

Write the spatial object as a shapefile with the updated Bicycle theft counts

writeOGR(boroughs, ".", "updated_boroughs", driver="ESRI Shapefile")

• Spatial autocorrelation
• Point pattern analysis
• Geographically weighted regression
• Kriging

Have a look at Bivand, Pebesma, & Rubio (2013) for further information.

Slides are available at
http://rpubs.com/pracademic/GIS_User_Forum