Mapping And Modelling Geographic Data In R

Session 3: Using R as a GIS

Introduction

The aim of this session is to give a flavour of how R can be used to store, manipulate and to analyse geographical data (data for which locations are recorded). For it, we shall focus on vector data, drawing extensively on the package documentation for Simple Features for R.

If you are interested in raster data, please take a look at the relevant parts of Geocomputation with R and the introduction to the raster package at rspatial.org.

Why use R as a GIS?

There are many other GIS software, including ArcGIS and the free, cross-platform and open-source QGIS. Why complicate matters by using R, which introduces the need to code and hasn’t the sorts of drop-down menus and other graphical user interfaces found in those other packages?

The answer, for me, concerns integration – the integration of a coding environment with a statistical environment, a graphical environment and so forth. It allows me to write scripts that process the data through all the stages of ‘data wrangling’ from inputting and collating the data, through analysing and fitting models, to communicating and mapping the results. That is why, for me, R provides an ideal environment for geographic data science and analysis.

In fact, it has been years since I last used a ‘proper’ GIS package because everything I need is in one of the many spatial packages for R, many of which are gathered here:

browseURL("https://cran.r-project.org/view=Spatial")

There are lots of them and their number is increasing! It is worth noting that there is a short cut way to installing many of these packages all at once. I have included the code to do so below but with the commands ‘commented out’. There is no need for you run the code here but it can be useful if you wish to install the packages on your own computer in the future.

## Do Not Run
# install.packages("ctv")
# ctv::update.views("Spatial")

Other packages are also grouped by Task View:

browseURL("https://cran.r-project.org/web/views/")

From sp to sf

Up until recently, the sp package (short for spatial) was the foundation for handling spatial data in R, providing classes and methods supporting geographic data, primarily vector data. This is being superseded by sf (short for simple features) that better meets the formal standard (ISO 19125-1:2004) for “how objects in the real world can be represented in computers, with emphasis on the spatial geometry of these objects” (https://cran.r-project.org/web/packages/sf/vignettes/sf1.html). At the time of writing, both are in operation and likely to be so for some time because of how many other packages build upon sp. There will be elements of both in this unit as there are also in books such as Brunsdon’s and Comber’s An Introduction to R for Spatial Analysis & Mapping (2nd edition).

What are vector data?

Vector data are those that are reducible to points. They are usually encoded with a \((x, y)\) coordinate but it could also include height \((x, y, z)\) and even a fourth dimension such as measurement accuracy.

if(!"sf" %in% installed.packages()[,1]) install.packages("sf")
if(!"tmap" %in% installed.packages()[,1]) install.packages("tmap")
require(sf)
require(tmap)
pt <- st_point(c(-2.5879, 51.4545))
pt <- st_sfc(pt, crs = 4326)
tmap_mode("view")
tm_basemap("Stamen.Watercolor") +
  tm_shape(pt) +
  tm_tiles("Stamen.TonerLabels") +
  tm_bubbles(col = "red")

Note that a preview of the various tm_basemaps is available at:

browseURL("http://leaflet-extras.github.io/leaflet-providers/preview/")

In any case, and regardless of the choice of basemap, the object, pt, is of class simple features column (sfc) and contains the geometry of the object (here just a single point feature) and also information about the coordinate reference system (crs):

pt

The coordinate reference was set and can be identified by its EPSG code – you can look it up using the website below:

browseURL("https://www.epsg-registry.org/")

Use the website to find the EPSG code for the British National Grid system.

A line is simply a collection of points, as in this .gpx file of a run route in Lyde Green (to the North East of Bristol):

if(!"tmaptools" %in% installed.packages()[,1]) install.packages("tmaptools")
require(tmaptools)
run_route <- read_GPX("Figure_Of_Height.gpx", layers = "tracks")
run_route
tmap_mode("view")
tm_basemap("OpenTopoMap") +
  tm_shape(run_route) +
  tm_lines(col = "red", lwd = 4)

A polygon is a collection of points where the start point is the same as the end point. As with points and lines, you can create them yourself…

pts <- rbind(c(1, 5), c(2, 2), c(4, 1), c(4, 4), c(1, 5))
multipoints <- st_sfc(st_multipoint(pts))
multipoints
a_line <- st_sfc(st_multilinestring(list(pts)))
a_line
a_polygon <- st_sfc(st_polygon(list(pts)))
a_polygon

… but it is more likely that you will read them in from other sources:

postcode_sectors <- st_read("england_pcs_2012.shp")
postcode_sectors
tmap_mode("plot")
  tm_shape(postcode_sectors) +
  tm_polygons(col = "red")+
  tm_compass(position = c("right", "top")) +
  tm_scale_bar(position = c("left", "bottom")) +
  tm_layout(title = "BS Postcode sectors", title.size = 0.85) +
  tm_logo("logo.png", position = c("right", "bottom"), height = 2,
          margin = 0.01)

Although tmap is often convenient, we don’t have to use it to plot an sf object. The plot command can be used although it doesn’t necessarily behave as we might want it to:

plot(postcode_sectors)

We can infer what is happening by looking at the data

if(!"tidyverse" %in% installed.packages()[,1]) install.packages("tidyverse")
require(tidyverse)
postcode_sectors %>%
  st_drop_geometry(.) %>%
  glimpse(.)

To just map the boundaries – the geometry – of the map…

postcode_sectors %>%
  st_geometry(.) %>%
  plot(.)

… which can then be modified in various ways. For example,

postcode_sectors %>%
  st_geometry(.) %>%
  plot(., col = sf.colors(12, categorical = TRUE), border = "white")
plot(st_geometry(st_centroid(postcode_sectors)), add = TRUE, pch = 21,
     bg = "white", cex = 0.5)

Try modifying the code above so that the centroids are represented by dark red squares. You may find the help menu ?points useful to achieve this – scroll down until you can see the ‘pch’ values.

We could also use ggplot2 for mapping. For example,

if(!"ggplot2" %in% installed.packages()[,1]) install.packages("ggplot2")
require(ggplot2)
centroids <- st_centroid(postcode_sectors)
ggplot() +
  geom_sf(data = postcode_sectors, fill = "red", col = "white") +
  geom_sf(data = centroids)

Another package to be aware of is ggmap as this can be useful for mapping sf objects over a Google basemap. In principle, the following code chunk will produce a .png file showing the run route from earlier layered on top of a Google road map.

if(!"ggmap" %in% installed.packages()[,1]) install.packages("ggmap")
require(ggmap)
png("road_map.png")
basemap <- get_map("bristol", maptype = "roadmap", source = "google",
                   zoom = 11)
plot(st_transform(st_geometry(run_route), crs = 3857)[1], bgMap = basemap,
     expandBB = rep(1.5, 4), col = "red", lwd = 2)
dev.off()

The output should look like:

browseURL("road_map.png")

However, to obtain the Google basemap does require that you have a registered API key. You can obtain this by signing-up at,

browseURL("https://cloud.google.com/maps-platform/")

– click on Get Started. Once you have the key, you can register it,

# register_google(key = "[Your key goes here]", write = TRUE)
## NB I have 'commented out' the code above to prevent it from being run
## You would run it without the initial # and with you key replacing
## [Your key goes here]

If you have registered your API then modify the code above to save the map as a jpeg file instead of a png, and have the run route later on top of a terrain mean. Try the help menus, ?png and ?get_map if you are not sure how to do it.

There is a very useful introduction to plotting simple features using the various approaches at https://r-spatial.github.io/sf/articles/sf5.html.

Simple Features for R

For the rest of this session primarily we will work through material developed by Edzer Pebesma as part of the simple features documentation. To begin, however, read through Chapter 2 of Lovelace et al. up to but not including Section 2.3:

browseURL("https://geocompr.robinlovelace.net/spatial-class.html")

Manipulating Simple Feature Geometries

The following code chunk opens the raw markdown file written by Edzer – it includes all the various blocks of code contained in the documentation. If you prefer to use the web document, it is at https://r-spatial.github.io/sf/articles/sf3.html but it will require cutting and pasting between it and the R Console.

download.file("https://raw.githubusercontent.com/r-spatial/sf/master/vignettes/sf3.Rmd", "sf3.Rmd")
file.edit("sf3.Rmd")

Manipulating Simple Features

https://r-spatial.github.io/sf/articles/sf4.html

download.file("https://raw.githubusercontent.com/r-spatial/sf/master/vignettes/sf4.Rmd", "sf4.Rmd")
file.edit("sf4.Rmd")

If you still have time, then work through…

Reading, Writing and Converting Simple Features

https://r-spatial.github.io/sf/articles/sf2.html

download.file("https://raw.githubusercontent.com/r-spatial/sf/master/vignettes/sf2.Rmd", "sf2.Rmd")
file.edit("sf2.Rmd")

Conclusion

This session has shown how R is able to adopt the same sort of functionality found in typical GIS software, with these most recently implemented through Simple Features for R. Before finishing, it is worth casting an eye over the list of functions contained in the online reference for sf. There also is a useful cheatsheet available.

Further Reading

Lovelace R, Nowosad J & Muenchow J (2019) Geocomputation with R, Chapters 2 to 5

Pebesma P (2018) Simple Features for R: Standardized Support for Spatial Vector Data PDF download. The R Journal 10 (1), pp. 439-46. https://journal.r-project.org/archive/2018/RJ-2018-009/RJ-2018-009.pdf