Introduction

In this handout we will see how to load the data for the spatial task of the assignment into R and visualise it using dynamic maps. These maps can’t be printed directly (although the output can be captured as screenshots). However they allow you to explore the data very thoroughly.

Starting the data analysis

Just as before the first step is to create a new project and a folder on the server to place the data.

Next you need to upload the zipped data into the folder. You will first have to download it from Brightspace (as you did when using QGIS). There is no need to unzip it this time, as the files will be automatically unzipped when they are uploaded to the server. So once you have the files ready go to the upload button in the pane on the bottom right of the Rstudio server and upload the files. You should now see them in the files pane.

We will carry out the analysis in several steps. The first step is to visualise the data, just as in any other type of analysis. However in this case the data will be visualised in the form of dynamic maps, which are rather similar to the maps you can produce in QGIS.

Although there are more options for exploring the data in depth in QGIS, one advantage of the maps produced by R is that they are standalone web pages. In other words there is no need for any software (apart from a web browser) in order to use them. Maps built in R can be sent to others as web pages or embedded in a webpage and will “work” just as they do on the server.

Making a markdown document

Once you have the project set up and the data loaded you need to start a new markdown document. You should be used to doing this now. Its really worth making sure that you alter the first chunk to switch off messaging as spatial analysis can throw up a lot that you don’t really want to include in your final document.

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

Libraries needed

The tidyverse packages are required for most work with spatial data. In order to load and analyse vector data the sf package is required. For some analyses (not those here) you also need the sp package, which uses the earlier format for spatial data. The mapview package is a really useful add on to R that allows spatial data to be viewed as “slippy” maps over online base maps. These maps can be zoomed, panned and interrogated using the mouse. Mapview maps are very similar to the interface provided to visualise data in a desktop GIS. The leafletextras package allows some extra features such as full screen controls and measuring tools to be added. The tmap package produces printable maps. Rcolorbrewer is useful for designing palettes and colour schemes for maps.

## This code chunk can be pasted at the start of any markdown document to add most of the key packages that makes R a GIS

library(tidyverse)
library(mapview)
library(sf)
library(tmap)
library(RColorBrewer)
library(elevatr)
library(raster)
library(units)
library(leaflet.extras)

Loading the data

The data are provided in the form of ESRI shapefiles. These are commonly used for spatial data due to the dominance of ESRI’s software for GIS work. Each data set is held as a combination of files with the same name, but with different extensions. You need to have all these files in order to work with the data (not just the file with .shp extension). However you only need to load in the file with this .shp extension in order to get the data into R.

birds<- st_read("Bird-data-v1.shp")
## Reading layer `Bird-data-v1' from data source 
##   `/home/rstudio/webpages/books/Bird_mapping_assignment/Bird-data-v1.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 16732 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 378500 ymin: 75500.02 xmax: 435500 ymax: 113000.1
## Projected CRS: OSGB 1936 / British National Grid

Notice the information provided by R about the data. They are points and the coordinate reference system used in British National grid.

habitat<-st_read("Land-cover-data-v1.shp")
## Reading layer `Land-cover-data-v1' from data source 
##   `/home/rstudio/webpages/books/Bird_mapping_assignment/Land-cover-data-v1.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 9 features and 1 field
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 378387 ymin: 75197.4 xmax: 435753 ymax: 113266
## Projected CRS: OSGB 1936 / British National Grid

In this case the geometries are multipolygons. In other words, there are multiple polygons associated with each land use type. It is possible to split these into single polygons for some analyses.

hwl<-st_read("High waterline.shp")
## Reading layer `High waterline' from data source 
##   `/home/rstudio/webpages/books/Bird_mapping_assignment/High waterline.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1065 features and 5 fields
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: 375137.1 ymin: 75196.2 xmax: 453490.5 ymax: 115446.4
## Projected CRS: OSGB 1936 / British National Grid

Mapview to look at the spatial data

The mapview package in R is a great tool for quickly looking at spatial data and planning analyses in R. When simply exploring data using mapview you are not necessarily interested in making a shareable map. Just as with other data a lot of visualisations are just part of the analytical work flow. When using spatial data we often just need to check that everything is suitable for the GIS operation that will produce some potentially useful numerical output.

However, mapview can also makes maps that are useful for sharing with others, as the maps take the form of dynamic webmaps in webpage

You do need to be aware of the limitations. Some things that can be done with a bespoke GIS such as QGIS are not possible using mapview.

  1. It cannot be used directly for very large, complex, spatial layers. As all the data has to be embedded in a webpage large layers will lead to large HTMLs. This may be OK up to a point but eventually they become too big to work as web pages. So look at the size and complexity of the layer before starting. If you are working with big data then use QGIS to visualise it, or for very advanced use R can be combined with a geoserver.
  2. It can be very tricky to find the right colours to visualise the data. R works with a wide range of palletes but while it is possible to specify individual colours for each category this is a time consuming process. Its easier to do with a point and click interface such as QGIS.
  3. Users of the web map may expect to be bale to interrogate the map in more detail than the interface provided is capable of doing. While extra elements (such as measurement tools) can be added users can expect more.
  4. Making a large number of individual mapviews in one markdown document can lead to confusion. A preliminary “exploratory” markdown can be useful for you to decide on how best to show the data, but then a simpler single map can be more useful for communicating with others.

Using mapview with a markdown document

There are a few things to be aware of when compiling markdown documents with mapview code chunks.

  1. Every map included in the final document adds to the size of the document, both in terms of space and in terms of the file size and speed at which the document can be rendered as a web page. It is not a good idea to include a lareg number of complex maps in the same document..
  2. Maps can be built using mapview using the same logic as ggplots. I.e. you can form a map as an object in R and add to it as you go on. You don’t need to make it visible until the end of the process.

Simplest possible mapview

Any one of the vector layers can be visualised by simply using it as an input to mapview. So two lines of R code gives a simple webmap (load the data then look at it in mapview). Here the line needed is commented out to prevent the handout becoming loaded down with too many maps. You can run it to see what it does then comment it out before knitting the document. …. Or run it in the console.

# mapview(habitat)

Setting the base maps

When the data are visualised using mapview a series of base maps are available. If you don’t set any options you will get default basemaps, but it is often useful to control this and select the base maps you want to see in the order you want to see them. You set this up before making the maps. Check for other options to use.

mapviewOptions(basemaps =c("Esri.WorldImagery",
                             "OpenStreetMap",
                              "Esri.WorldShadedRelief"), homebutton = FALSE)

Setting up a colour pallete.

Finding the best colours to use for any map can be the hardest part of the process. The RColorBrewer package has a range of options. Check the documentation to see them.

pal<- brewer.pal(12,"Paired")

Maps

Birds

The point data on the bird observations can be visualised by species.

If the “burst” option is used when forming a mapview then each species can be selected in turn. It is always worth adding a full screen option to the maps (from the leaflet.extras package) in order to look at the map more closely.

The palette set up above provides the colours. You can try other palettes

Notice that the map can be saved as an object in R before visualising it by typing its name. This can be useful for adding maps together.

mapview(birds,zcol ="Species",burst=TRUE, hide=TRUE, col.regions=pal) -> birdmap
birdmap@map %>% addFullscreenControl() -> birdmap@map

Again I have commented out the line that produces the map here to prevent the document becoming too heavy when compiled.

# birdmap

Land cover

The same approach can be used to produce a map of the land covers.

Finding the area of each land use type

An area attribute can be added to the data. In this case the area would be the total for all the polygons making up each land use class. The results are provided in the form of a “units” class, so the units library can be used for easy conversion between square meters and square kilometers. Note that as the geometries are very complex they need to be dropped before looking at the results as a table.

habitat$Area<-drop_units(set_units(st_area(habitat), km^2))

habitat %>% st_drop_geometry()
##    LandCover      Area
## 1     Arable 316.56388
## 2    Coastal  31.48573
## 3 Freshwater  11.08190
## 4  Grassland 479.09516
## 5  Heathland 170.76895
## 6  Saltwater  11.86812
## 7   Suburban 183.56570
## 8      Urban  30.44617
## 9   Woodland 330.05616
mapview(habitat,zcol="LandCover",burst=TRUE,hide=TRUE, col.region=pal) -> habitatmap
habitatmap@map %>% addFullscreenControl() -> habitatmap@map
# habitatmap

Combined

The two previous maps can be added together to produce a combined map.

It can be useful to add a measuring tool to maps. This is also from the leaflet extras package.

Notice that there is what seems to be a bug in the mapview package here. When the maps are added all the options for the bird species appear as legends even though only one is shown. If you click through all the species in turn this is resolved.

birdmap + habitatmap -> combinedmap
combinedmap@map %>%  addMeasure(primaryLengthUnit ='meters', secondaryLengthUnit='kilometers',primaryAreaUnit='hectares',position="topleft")