Chapter 3 Some simple mapping

3.1 Introduction

This is a very quick demonstration of how R can be used to make maps that are embedded in markdown documents. The data used is stored on the server in this case, so the code will not work locally. However the data layers could be saved as shapefiles and exported in the usual way for local use off the server.

3.2 Getting started

1. Open a new project. You might want to call it Hengistbury as this is the area which will be the focus of mapping.
2. Start a new markdown document and save it with an appropriate name. Delete all the material below the initial chunk. You can very carefully modify the first chunk to prevent annoying messages being added to your document.

knitr::opts_chunk$set(echo = TRUE, warning = FALSE,message = FALSE)

3.3 Load some libraries

The tidyverse library should always be loaded at the start of any analysis, as without this R cannot run many common functions that are used for “tidy” data handling. To turn R into a GIS requires additional libraries to be loaded.

library(tidyverse)
library(giscourse)
library(mapview)
mapviewOptions(basemaps =c("Esri.WorldImagery","OpenStreetMap", "CartoDB.Positron", "OpenTopoMap"))
library(tmap)
tmap_options(check.and.fix = TRUE)
library(sf)

3.4 Load the data

The data for this exercise is contained within the giscourse package that was loaded in the previous chunk.

data(hengistbury)

3.5 Mapview

The “lnr” object that can be seen in the Environment tab is the limits of the local nature reserve. This can be turned into a slippy web map with the following code.

mapview(lnr,  alpha.regions = 0, legend=FALSE, lwd =3) %>% extras()

3.6 Classic paper map

A shaded “chloropleth” map of the priority habits maps provided by natural England can be made using the tmap package. This uses styles and themes to decide on the colours used. Finding the correct colours for a map can be a time consuming process. The map was sourced from here before loading into the server. https://data.gov.uk/dataset/4b6ddab7-6c0f-4407-946e-d6499f19fcde/priority-habitat-inventory-england

habmap<-qtm(habitats, "mhabitat") +tm_style("classic") + tm_scale_bar(position="left")
habmap

This map can also be viewed in “interactive” mode.

tmap_mode("view")
habmap

The following code produces a mapview in which each habitat shows as a separate layer.

mapview(habitats, zcol="mhabitat", burst=TRUE,legend=FALSE,  hide = TRUE) %>% extras()

3.7 Simple digitising

On screen digitising is best conducted using a dedicated desktop GIS such as QGIS. However, if you just want to quickly draw around some features or place some points on a map it is easy to do directly in R. Digitising cannot be done within a markdown document, so the function must be run in the console.

Type “digitise()” into the console. You will see a map appear in the plots window. Zoom this to full screen and use the tool to draw either points or polygons. I very quickly drew around some features representing heathland and grassland in the first field within the reserve. As a demonstration I just used four polygons and did not take time to draw around the shapes very carefully. It is easy to do a much better job, even in R! However this shows the concept.

mymap<-st_read("mymap.shp")

## Reading layer `mymap' from data source 
##   `/home/rstudio/webpages/books/qs_21/mymap.shp' using driver `ESRI Shapefile'
## Simple feature collection with 4 features and 3 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -1.768694 ymin: 50.7164 xmax: -1.763546 ymax: 50.71906
## Geodetic CRS:  WGS 84

mapview(mymap,zcol= "id", burst=TRUE)

You should notice that the polygons have been drawn over each other. If we want to produce a layer representing different habitats this is not correct, as it suggests that an area on the map can be both grassland and heathland at the same time. It is quite simple to fix this in R using a differencing operation.

mymap<-st_difference(mymap)
mapview(mymap,zcol= "id", burst=TRUE)

3.8 Adding attributes

When digitising in a desktop GIS the attributes are added during the drawing phase. This is not possible in R. However each polygon has a unique id which can be merged with a data frame containing attributes. This data frame can be a csv file that is produced in Excel. So if you know the id number of a polygon or point on the map and have data that corresponds to that feature you can add it to the shapefile.

id<-1:4
habitat<-c("heath", "heath", "heath", "grassland")
d<-data.frame(id,habitat)
mymap<-merge(mymap,d)

3.9 Calculating areas

mymap$area<-st_area(mymap)

mapview(mymap,zcol= "habitat", col.regions=c("darkgreen", "brown"))

3.10 More on digitising

Although the digitising capabilities of R are rather more limited than a desktop GIS, the sf package in R has many powerful functions built in that can run basic geoprocessing operations. For example let’s digitise some paths at Hengistbury head as lines.

# Run in console
digitise(shp="paths.shp")

Now we can read in the shapefile we produced.

paths<-read_sf("paths.shp")

The map has been digitised as latitude and longitude, as can been seen by checking the coordinate reference system.

st_crs(paths)

## Coordinate Reference System:
##   User input: WGS 84 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["latitude",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["longitude",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4326]]

This can be changed easily in R to British national grid.

paths<-st_transform(paths,27700)

Now we can try some simple operations such as buffering the lines to find the area within ten meters of a path.

paths_buffer<-st_buffer(st_union(paths), 10)
mapview(paths_buffer)