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 datavery thoroughly.

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.

library(tidyverse)
library(mapview)
library(sf)
library(tmap)
library(RColorBrewer)
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/qs24_spatial/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/qs24_spatial/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.

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.

library(units)
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

Setting the base maps

When the data are visualised using mapview a series of base maps are available. The basemaps can be set 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. 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()

Land cover

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

mapview(habitat,zcol="LandCover",burst=TRUE,hide=TRUE, col.region=pal) -> habitatmap
habitatmap@map %>% addFullscreenControl()

Combined

The two previous maps can be added together to produce a combined map. 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.

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

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

Saving the workspace

The workspace, which now contains both the data and the maps that have been built from it can be saved to a file. This can be useful when several markdown documents are used. The work done in the first document of the series can then be passed onto the next one. I tend to use the same name as the document for the saved workspace, as this makes it clear where the workspace came from. However, be very careful not to accidentally use an Rmd extension for the workspace, as this would overwrite the markdown file itself. Use .rda for the saved data.

save(list=ls(),file="Exploratory_maps.rda")