Although the zoomable maps produced by mapview provide the most powerful tools for investigating the data yourself, there will often come a point when you are writing up an analysis when you need to produce maps as figure that can be placed in the article.
It can be very challenging to show complex spatial information in a clear and communicable manner. Maps tend to become very “busy” and it is quite easy to lose their context. In many cases it is almost impossible for the reader to extract much useful information from visual inspection of a printed map in the form of an analysis. In reality this often does not matter too much, as it is the job of the GIS to do things such as overlaying the points on the polygons and extracting the information, not the person looking at the map.
The tmap package is probably the best way to make printable maps although there are many other options.
library(tidyverse)
library(mapview)
library(sf)
library(tmap)
library(RColorBrewer)
library(rnaturalearth)
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
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
This is a very quick “default” tmap. The only complication here is the filter used to extract a single bird species from the data.
Note how elements such as scale bar and compass are added to the map.
You can explore the other themes available to see how the look of the map is altered. T
tmap_style("albatross")
qtm(habitat, fill="LandCover") +qtm(filter(birds,Species=="Yellowhammer"), title="Yellowhammer", fill="red") + tm_scale_bar() +tm_compass(type="4star") ->yh
yh
A map that helps the reader to locate the study site can be useful. This can be placed in the methods section, rather than the results.
The rnaturalearth package holds some rather low resolution outline maps of countries and administrative regions (in the case of the UK these are counties). These can be useful as a simple base when plotting the spatial data.
We can get the UK counties with this line of code. Notice the use of st_transform. The data are provided with coordinates as latitude and longitude. But we are working with data in the British National grid coordinate reference system. This is EPSG code 27700. So the data are reprojected into this CRS.
ne_states(country="United Kingdom", return="sf") %>% st_transform(crs=27700)-> UK
We might just want a subset of counties. To get this you can use a filter that looks for matches in a lost of names. The names have to be precisely the same as those in the data.
UK %>% filter(name %in% c("Dorset","Hampshire", "Isle of Wight")) -> South
To help the reader locate the study are we might start with a map showing where the southern counties lie in the UK.
tmap_style("natural")
qtm(UK, style ="natural", title="UK") +qtm(South,fill="red") ->ukmap
ukmap
We could then add the habitat (land use) map to just the southern counties.
tmap_style("classic")
qtm(South) +qtm(habitat,legend=FALSE, title="Land use") ->hmap
hmap
To give an impression of the number of point observations of all bird species we could combine the counties with the bird data.
qtm(South) +qtm(birds,legend=FALSE, title="bird observations") ->bmap
bmap
Or do this just for one species. Note the use a filter to select the species.
qtm(South) +qtm(filter(birds,Species=="Yellowhammer"),legend=FALSE, title="Yellowhammer observations") ->yhmap
The four maps produced by the code above can be placed together as a single figure using the tm_arrange function.
tmap_arrange(ukmap,hmap,bmap,yhmap, asp=NA, ncol=2) -> fig1
fig1