Mapping with QGIS and R

Justin Dyck, MSc

2022-10-24

QGIS

  • QGIS is a professional GIS application that is built on top of and proud to be itself Free and Open Source Software (FOSS).
  • This is an alternative to closed source software like ESRI’s ArcGIS which costs $845/year to start.
  • QGIS does the job needed for this class, and can do a lot more if GeoStatistical tools are needed for future research. I will only cover mapping here, but there is a lot of functions that go way beyond this.

https://www.qgis.org/en/site/

Starting with QGIS

  • QGIS will work with many different file types.
  • For this purpose, we will be using Shape files (.shp extension), and .csv data files.
  • The Shape file creates the polygonal map, and the .csv holds the output data (from whatever model output we have).
  • We’ll go through loading the Shape file, linking the data, and then displaying it to print to an image file.

Loading Data

  • First the Shape file:
  • Layer -> Add Layer -> Add Vector Layer
  • Choose the File path and file in the “Source” section

Loading Data

  • This gets us an empty map:

Loading Data

  • Now we need to load our .csv data
  • Layer -> Add Layer -> Add Delimited Text Layer
  • Choose the file from the File name option
  • Ensure to check the “No geometry” option in the Geometry Definition section

Linking Data

  • To join the 2 layers, right click on the Map layer in the “Layers” tray (bottom left of screen). Choose the “Properties” option.
  • This brings up all the options for this layer.
  • Navigate to the Joins page.
  • Click on the ‘+’ symbol.
  • Now choose the layer to join to the map, and the fields you want them to join on. (Here would be the “Number” fields).

Displaying Data

  • Open up Map Layer’s properties again
  • Navigate to the “Symbology” page.
  • At the top select “Graduated” then select the variable to map in the “Value” option

Displaying Data

  • We then select the required mode, and hit “Classify”.
  • Select different palette’s from the “Color Ramp” option
  • Then “Apply”
  • Can also choose different classes, and other options to tweak and customize the map aesthetics

Displaying Data

Displaying Data

Mapping with R

  • R has many spatial and GIS packages currently available.
  • Also, many ways to display data.
  • We will focus on the sf package for data manipulation, and the tmap and leaflet packages for displaying maps.
library(sf)
library(tmap)
library(leaflet)

sf_use_s2(FALSE)

Loading Data into R

  • While there are ways (using rgdal is one way), I will use the sf built in functions.
  • This provides an easy way to manipulate the spatial data.
  • For later, we also need to set the coordinate system. This will allow us to overlay the shapefile onto larger basemaps.
MB_Shape = st_read(("Manitoba_map67.shp")) %>%  #Read shapefile
  st_transform(crs = 4326) #Transform coordinate system
Reading layer `Manitoba_map67' from data source 
  `/media/justin/Extension/Work Projects/Data Science/Workshops/SpatialEpi_Lecture/Manitoba_map67.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 67 features and 4 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -102.0123 ymin: 48.99876 xmax: -88.98969 ymax: 60.00171
Geodetic CRS:  GCS_Assumed_Geographic_1

Manipulating Spatial Data

head(MB_Shape)
Simple feature collection with 6 features and 4 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -97.27762 ymin: 48.99875 xmax: -95.15352 ymax: 50.28392
Geodetic CRS:  WGS 84
          FIRST_NAME           RHA       NAME2 Number
1  BN2-20 Brokenhead North Eastman  Brokenhead     11
2 BN5-20 Springfield North Eastman Springfield     13
3     BS1-25 Central South Eastman     Central     16
4    BS2-25 Northern South Eastman    Northern     17
5  BS3-25 South/East South Eastman    Southern     18
6     BS4-25 Western South Eastman     Western     19
                        geometry
1 MULTIPOLYGON (((-96.35287 5...
2 MULTIPOLYGON (((-96.70431 5...
3 MULTIPOLYGON (((-96.64578 4...
4 MULTIPOLYGON (((-96.36583 4...
5 MULTIPOLYGON (((-97.1671 49...
6 MULTIPOLYGON (((-96.98073 4...

Manipulating Spatial Data

  • We also need to add in data (what we want to map)
RR = read_csv(("ExampleData.csv")) %>% 
  mutate(Number = as.character(Number)) #This is an index not a numeric variable

JoinedData = left_join(

  MB_Shape %>% mutate(Number = as.character(Number)),
  RR,
  
  by = c("Number" = "Number")
)

Mapping with tmap

tm_shape(JoinedData) +  #tm_shape starts the map
  tm_polygons("RR",     #tm_polygons tells it to build a map using the polygons in the shapefile and fill them with our data
              palette = sequential_hcl(12, "reds 3", rev = T),
              title = "Risk Ratio")

Mapping with tmap

  • Can also use different breaks:
tm_shape(JoinedData) + 
  tm_polygons("RR",     
              palette = sequential_hcl(12, "reds 3", rev = T),
              title = "Risk Ratio",
              breaks = c(quantile(JoinedData$RR, probs = seq(0, 1, by = 0.2))))

Mapping with tmap

  • And many other options!
Map = tm_shape(JoinedData) + 
  tm_polygons("RR",     
              palette = sequential_hcl(12, "reds 3", rev = T),
              title = "Risk Ratio",
              breaks = c(quantile(JoinedData$RR, probs = seq(0, 1, by = 0.2)))) +
  tm_scale_bar(position = c("left","BOTTOM"))+      #Add a scale bar
  tm_layout(title = "Example in MB",            #layout options to make things pretty and annotate if necessary
            legend.position = c("right","bottom"),
            inner.margins = c(0.05,0.05,0.12,0.05),
            frame=T,
            legend.show=T,
            bg.color = "white",
            legend.bg.color = "white",
            legend.frame = T)

Mapping with tmap

Multiple Maps

  • Maybe we want to zoom and see Winnipeg more clearly.
  • Maybe we want to see both MB and Winnipeg maps side by side!
WPG_Data = JoinedData %>% 
  filter(RHA == "Winnipeg")

WPG_Map = tm_shape(WPG_Data) + 
  tm_polygons("RR",     
              palette = sequential_hcl(12, "reds 3", rev = T),
              title = "Risk Ratio",
              breaks = c(quantile(JoinedData$RR, probs = seq(0, 1, by = 0.2)))) +
  tm_scale_bar(position = c("left","BOTTOM"))+     
  tm_layout(title = "Example in WPG",           
            legend.position = c("right","bottom"),
            inner.margins = c(0.05,0.05,0.12,0.05),
            frame=F,
            legend.show=F, #Hide legend here
            bg.color = "white",
            legend.bg.color = "white",
            legend.frame = T)

Multiple Maps

tmap_arrange(Map, WPG_Map)

Interactive Maps With Leaflet

  • Leaflet is a powerful JavaScript mapping library
  • Used in many web-based mapping applications
  • R has a leaflet package, that utilizes this JS library to allow us to have the same functionality within R.

Interactive Maps With Leaflet

  • First we create a basemap and set the zoom:
leaflet() %>% 
  addProviderTiles(providers$CartoDB.Positron) %>% # background map
  setView(lng =  -97.8, lat = 54.6, zoom = 5)

Interactive Maps With Leaflet

  • Next, we create a palette, and add in the data
pal = colorNumeric(
  palette = sequential_hcl(12, "reds 3", rev=T),
  domain = JoinedData$RR
) 

LeafletMap = leaflet() %>% 
  addProviderTiles(providers$CartoDB.Positron) %>%
  setView(lng =  -97.8, lat = 54.6, zoom = 5) %>% 
  addPolygons(data = JoinedData, weight=1, color="#000000",
              opacity=0.4,
              fillOpacity = 0.6, smoothFactor = 0.1,
              fillColor = ~pal(JoinedData$RR)
  ) %>%
  addLegend("bottomleft", pal = pal, values = quantile(JoinedData$RR,na.rm=T),
            title = "Risk Ratio",
            opacity = 1
  )

Interactive Maps With Leaflet

Publishing

  • Not the scope of this workshop, but publishing of any maps can be done within R easily!

- Saving a tmap with the jpeg()/tiff() functions, in the same way you would save a plot image to paste into a manuscript or report.
- Using rmarkdown to write a pdf manuscript or report. tmap’s can be called in line (similar to calling an image in latex, but using R-code rather than latex).
- html reports/vignettes/slides supports both tmap or interactive html images like leaflet!
- The flexdashboard package (https://rstudio.github.io/flexdashboard/articles/using.html), is a great way to use rmarkdown to create simple dashboards. This will support static or interactive maps.
- The Shiny package (https://shiny.rstudio.com/). This is a powerful tool to create complex dashboards and interactive data applications. The learning curve is a bit higher here, but it is well worth the investment if looking for a job outside of academia.

Additional Resources

  • sf: https://r-spatial.github.io/sf/articles/
  • tmap: https://cran.r-project.org/web/packages/tmap/vignettes/tmap-getstarted.html
  • leaflet: https://rstudio.github.io/leaflet/