GIS in R (with NZ data)

Isaac Bain

2017-04-01

R is an ever increasingly powerful option for conducting GIS-type analyses and easily producing static and interactive maps. Reproducability is a key attribute for time saving and open science.

Here I provide some examples of where to obtain, and how to work with, New Zealand data. In this tutorial we will download the point locations of all 961 DOC huts and place them on an interactive, lightweight, webmap.

R packages

In this tutorial I will be using leaflet to produce interactive maps. These other packages are neccessary for loading and manipulation of spatial data.

library(tidyverse)
library(sp)
library(leaflet)
library(spdplyr)
library(rgdal)

Data sourcing

New Zealand is fortunate to have relatively free and open access to goverment spatial data. A good place to begin searching for spatial data is at Koordinates. Similar data services are maintained by the likes of Ministry for the Environment, LINZ, etc.

The most common contemporary projection for working with NZ data is NZGD2000 / New Zealand Transverse Mercator 2000 (NZTM) - EPSG:2193. It is a good habit to download and work with data in this projection where possible, especially if you might be collaborating with other GIS users.

DOC data

In this tutorial we will be working with the DOC huts layer, which simply shows points where all the DOC huts are located. For your convenience I have uploaded both these files to my website so they are directly accessible from R, but this would work the same if you downloaded from Koordinates yourself.

download.file("http://www.isaacbain.co.nz/doc-huts.zip", destfile = "doc-huts.zip")  # downloads layer to your working directory

unzip("doc-huts.zip")  # unzips file, searching in working directory

Reading and re-projecting spatial data

By far the easiest method for reading in spatial data is when the data is obtained in shapefile format. This means that your projection is already defined (in a way that R can read), and there is no need to tell R which columns contain x and y coordinates.

Leaflet requires all data to be in the WGS84 projection, in order to match the common basemaps it provides. Fortunately it is easy to reproject spatial data using spTransform; in our case from NZTM2000 to WGS84.

huts <- readOGR("doc-huts","doc-huts")  # first argument is directory, second agrument is layer name

wgs84 = '+proj=longlat +datum=WGS84'  # define wgs84 information
huts <- spTransform(huts, CRS(wgs84))  # reproject to wgs84

Making maps

If you are a current user of dplyr or ggplot2 then making maps with leaflet will be very familar to you. It combines elements of both these packages; with its use of the pipe operator, %>% and the grammar of graphics approach to producing graphics by incrementally adding layers to a plot.

Lets make our first plot, of all the DOC huts in NZ.

  1. Tell leaflet where our spatial data is located, we called it huts.
  2. Add red points to our map to represent the huts. As at certain scales many points overlap, make them half transparent. addMarkers are ugly so I avoid them.
  3. Define the view we are most interested in, if the default view is unacceptable.
  4. Add a pretty basemap.
leaflet(huts) %>%  # map of all huts
  addCircleMarkers(stroke = FALSE, radius = 2.5, color = "Red", fillOpacity = 0.5) %>%
  setView(lng = 173.28 , lat = -41.27, zoom = 5) %>%  # centered around Nelson
  addTiles()  # default tiles are Open Street Map

Adding popups

Popups are a useful addition to many types of maps, to allow the user to glean extra information from the attribute table. Here we will add a popup = argument to the points layer.

leaflet(huts) %>%  # map of all huts
  addCircleMarkers(stroke = FALSE, radius = 2.5, color = "Red", fillOpacity = 0.5, 
                   popup = huts$DESCRIPTIO) %>%
  setView(lng = 173.28 , lat = -41.27, zoom = 5) %>%  # centered around Nelson
  addTiles()  # default tiles are Open Street Map

We can also add a few lines of information to popups and format them nicely using basic html tags. Here we first define how we want the popups to look, and then refer to this in the popup argument. At the same time we have added a label argument which allows users to rapidly view information without having to click.

state_popup <- paste0("<strong>Name: </strong>", huts$DESCRIPTIO, "<br><strong>Type: </strong>", huts$OBJECT_TYP, "<br><strong>Status: </strong>", huts$STATUS)

leaflet(huts) %>%  # map of all huts
  addCircleMarkers(stroke = FALSE, radius = 2.5, color = "Red", fillOpacity = 0.5, 
                   popup = state_popup, label = huts$DESCRIPTIO, labelOptions = labelOptions(textOnly = TRUE)) %>%
  setView(lng = 173.28 , lat = -41.27, zoom = 5) %>%  # centered around Nelson
  addTiles()  # default tiles are Open Street Map

Adding basemaps

Adding basemaps relevant to the data you wish to display is a neccessary component of cartography. Leaflet includes many useful default options, a complete list can be found here. The Stamen series is my favourite, and the Stamen Terrain is particulary nice.

It is also possible to add any WMS basemap to these maps, and luckily LINZ supplies many of there layers in this format. This has the big advantage of being able to serve large basemaps with the hosting being done by somebody else!

leaflet(huts) %>%  # map of all huts
  addCircleMarkers(stroke = FALSE, radius = 2.5, color = "Red", fillOpacity = 0.5, 
                   popup = state_popup, label = huts$DESCRIPTIO, 
                   labelOptions = labelOptions(textOnly = TRUE)) %>%
  setView(lng = 173.28 , lat = -41.27, zoom = 5) %>%  # centered around Nelson
  addTiles() %>%   # default tiles are Open Street Map
  addProviderTiles("Stamen.Toner", group = "Toner") %>%
  addProviderTiles(providers$Stamen.Terrain, group = "Terrain") %>% 
  addWMSTiles(
    "http://wms.data.linz.govt.nz/0f213e5b98d94b09b79f91f7cee570f7/r/wms",
    layers = "x767",
    options = WMSTileOptions(format = "image/png", transparent = TRUE),
    attribution = "LINZ CC BY 3.0",
    group = "LINZ 1:50 Topo") %>% 
  addLayersControl(
    baseGroups = c("Terrain", "Toner", "LINZ 1:50 Topo","OSM"),
    options = layersControlOptions(collapsed = FALSE))